##########################################################################################
#Title: Replication Code for "Psychology, Soft Skills, or Cash? Evidence on Marginal Investments"
#Author: Megan Lang
#Date prepraed: August 22, 2023
##########################################################################################

library(ggplot2)
library(readxl)
library(dplyr)
library(foreign)
library(gridExtra)
library(lfe)
library(car)
library(stringr)
library(stargazer)
library(xtable)
library(devtools)
library(extrafont) 
library(fwildclusterboot)

set.seed(20230618)
dqrng::dqset.seed(20230618)

##########################################################################################
#Read in the data
##########################################################################################
bl <- read.csv("bl.csv")
fu <- read.csv("fu.csv")
fu.lp <- read.csv("fulp.csv")
pw <- read.csv("pw.csv")
coop <- read.csv("Co-op.csv")
wksp <- read.csv("Wksp.csv")

##########################################################################################
#Construct Outcome Variables and Build Functions
##########################################################################################
#####Function for generating indices#####
gen.ind <- function(cols, col.names){
  bl.ind <- select(bl, resp_id, all_of(cols))
  fu.ind <- select(fu, resp_id, all_of(cols))
  ncol <- length(cols) + 1
  
  for (i in 1:length(bl.ind$resp_id)){
    raw <- as.matrix(bl.ind[i, c(2:ncol)])
    bl.ind$ewindex[i] <- mean(raw, na.rm=T)
  }
  
  for (i in 1:length(fu.ind$resp_id)){
    raw <- as.matrix(fu.ind[i, c(2:ncol)])
    fu.ind$ewindex[i] <- mean(raw, na.rm=T)
  }
  
  bl.ind <- select(bl.ind, resp_id, bl.indew = ewindex)
  fu.ind <- select(fu.ind, resp_id, fu.indew = ewindex)
  ind <- full_join(bl.ind, fu.ind, by = "resp_id")
  colnames(ind) <- c("resp_id",col.names)
  
  full <- left_join(full, ind, by = "resp_id")
}

#Function for post-workshop indices
gen.ind.pw <- function(cols, col.names){
  pw.ind <- select(pw, resp_id, all_of(cols))
  ncol <- length(cols) + 1
  
  for (i in 1:length(pw.ind$resp_id)){
    raw <- as.matrix(pw.ind[i, c(2:ncol)])
    pw.ind$ewindex[i] <- mean(raw, na.rm=T)
  }
  
  pw.ind <- select(pw.ind, resp_id, pw.indew = ewindex)
  colnames(pw.ind) <- c("resp_id",col.names)
  
  full <- left_join(full, pw.ind, by = "resp_id")
}

#####Function for performing multiple inference corrections####
# Note that this is taken directly from Michael Anderson's Stata code, which "generates BKY (2006) sharpened two-stage q-values 
# as described in Anderson (2008), "Multiple Inference and Gender Differences in the Effects of Early Intervention: A Reevaluation 
# of the Abecedarian, Perry Preschool, and Early Training Projects", Journal of the American Statistical Association, 103(484), 1481-1495

# The input to this function should be a dataframe of p-values from a given family, sorted from smallest to largest. The dataframe
# should have columns "Outcome", "coef", "pval", and "rank". Outcome is the name of the outcome variable. coef is the estimated
# coefficient on the treatment of interest for that outcome. pval is the naive p-value of the estimate.
# rank is the rownumber based on the sorted p-values.
gen.qvals <- function(df){
  totalpvalues <- as.numeric(length(df$Outcome))
  df <- mutate(df, qval = 1)
  #Set initial condition for the loop to calcualte q-values with 2-step procedure
  qv <- 1
  
  while (qv>0){
    #Generate ranking by p-value
    df <- arrange(df, pval) %>% mutate(rank = row_number())
    # First Stage
    # Generate the adjusted first stage q level we are testing: q' = q/1+q
    qval_adj <- qv/(1+qv)
    #Generate value q'*r/M
    df <- mutate(df, fdr_temp1 = qval_adj*rank/totalpvalues)
    #Generate binary variable checking condition p(r) <= q'*r/M
    df <- mutate(df, reject_temp1 = if_else(fdr_temp1>=pval,1,0))
    #Generate variable containing p-value ranks for all p-values that meet above condition
    df <- mutate(df, reject_rank1 = reject_temp1*rank)
    #Record the rank of the largest p-value that meets above condition
    total_rejected <- max(df$reject_rank1)
    #Second Stage
    #Generate the second stage q level that accounts for hypotheses rejected in first stage: q_2st = q'*(M/m0)
    qval_2st <- qval_adj*(totalpvalues/(totalpvalues-total_rejected))
    #Generate value q_2st*r/M
    df <- mutate(df, fdr_temp2 = qval_2st*rank/totalpvalues)
    #Generate binary variable checking condition p(r) <= q_2st*r/M
    df <- mutate(df, reject_temp2 = if_else(fdr_temp2>=pval,1,0))
    #Generate variable containing p-value ranks for all p-values that meet above condition
    df <- mutate(df, reject_rank2 = reject_temp2*rank)
    #Record the rank of the largest p-value that meets above condition
    total_rejected2 <- max(df$reject_rank2)
    #A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition
    df <- mutate(df, qval = if_else(rank <= total_rejected2, qv, qval))
    df <- select(df, Outcome, coef, pval, rank, qval)
    # Reduce q by 0.001 and repeat loop
    qv <- qv - .001
  }
  df <- mutate(df, adj.t_val = qt(qval/2, 453), 
               adj.SE = abs(coef/adj.t_val))
  return(df)
}

#####Construct pre-specified indices#####
#Create a dataframe with demographic information
demo <- select(bl, resp_id, workshop_name, age = d_24, ed = d_25, hhsize = d_26, men = d_27a,
               women = d_27b, girls = d_27c, boys = d_27d) 
#Self-Advocacy (do this one manually to check then use the function)
cols <- c("sei_12","sei_17a","sei_17b","sei_15")
bl.sa <- select(bl, resp_id, all_of(cols))
fu.sa <- select(fu, resp_id, all_of(cols))

for (i in 1:length(bl.sa$resp_id)){
  raw <- as.matrix(bl.sa[i, c(2:5)])
  bl.sa$ewindex[i] <- mean(raw, na.rm=T)
}

for (i in 1:length(fu.sa$resp_id)){
  raw <- as.matrix(fu.sa[i, 2:5])
  fu.sa$ewindex[i] <- mean(raw, na.rm=T)
}

bl.sa <- select(bl.sa, resp_id, bl.saew = ewindex)
fu.sa <- select(fu.sa, resp_id, fu.saew = ewindex)
sa <- full_join(bl.sa, fu.sa, by = "resp_id")

full <- left_join(demo, sa, by = "resp_id")

rm(demo, sa, bl.sa, fu.sa)

#Problem solving
ps.cols <- c("sei_5a","sei_5b")
ps.cn <- c("bl.psew","fu.psew")
ps.cn.pw <- c("pw.psew")
full <- gen.ind(ps.cols, ps.cn)
full <- gen.ind.pw(ps.cols, ps.cn.pw)

#Connections within Homes and Villages
chv.cols <- c("sei_10b", "sei_10d", "sei_11b", "sei_11d", "sei_11f")
chv.cn <- c("bl.conhvew","fu.conhvew")
full <- gen.ind(chv.cols, chv.cn)

#Connections within Cooperatives
cc.cols <- c("sei_10a","sei_10c","sei_11a","sei_11c","sei_11e")
cc.cn <- c("bl.conew","fu.conew")
full <- gen.ind(cc.cols, cc.cn)

#Communication skills
com.cols <- c("sei_6","sei_7a","sei_8a","sei_16a","sei_16b","sei_16c")
com.cn <- c("bl.comew","fu.comew")
full <- gen.ind(com.cols, com.cn)

#Self-Evaluation of Value
v.cols <- c("sei_3a","sei_3b")
v.cn <- c("bl.vew","fu.vew")
v.cn.pw <- c("pw.vew")
full <- gen.ind(v.cols, v.cn)
full <- gen.ind.pw(v.cols, v.cn.pw)

#Resilience
r.cols <- c("sei_4a","sei_4b")
r.cn <- c("bl.rew","fu.rew")
r.cn.pw <- c("pw.rew")
full <- gen.ind(r.cols, r.cn)
full <- gen.ind.pw(r.cols, r.cn.pw)

#####Leadership (# positions and holding any position)#####
bl.l <- select(bl, resp_id, sei_9__0, sei_9__1, sei_9__2, sei_9__3,
               sei_9__4, sei_9__5, sei_9__6, sei_9__7, sei_9__8, sei_9__9)
bl.l <- mutate(bl.l, nlead = if_else(sei_9__0=="NONE"|sei_9__0=="None"|sei_9__0=="Non"|
                                       sei_9__0=="ntayo"|sei_9__0=="Ntayo"|sei_9__0=="N/A"|
                                       sei_9__0=="oya",0,
                                     if_else(sei_9__1=="##N/A##",1,
                                             if_else(sei_9__2=="##N/A##",2,
                                                     if_else(sei_9__3=="##N/A##",3,
                                                             if_else(sei_9__4=="##N/A##",4,
                                                                     if_else(sei_9__5=="##N/A##",5,
                                                                             if_else(sei_9__6=="##N/A##",6,
                                                                                     if_else(sei_9__7=="##N/A##",7,
                                                                                             if_else(sei_9__8=="##N/A##",8,
                                                                                                     if_else(sei_9__9=="##N/A##",9,as.double(NA))))))))))),
               lead.ind = if_else(nlead>0,1,0))

fu.l <- select(fu.lp, resp_id = Respondent.s.ID.number, nlead = lpcount)
fu.l <- mutate(fu.l, lead.ind = if_else(nlead>0,1,0))

bl.l <- select(bl.l, resp_id, bl.nlead = nlead, bl.lead.ind = lead.ind)
fu.l <- select(fu.l, resp_id, fu.nlead = nlead, fu.lead.ind = lead.ind)
lead.ind <- full_join(bl.l, fu.l, by = "resp_id")

full <- left_join(full, lead.ind, by = "resp_id")

rm(bl.l, fu.l, lead.ind)

#####Self-Perceived Social Status####
bl.ss <- select(bl, resp_id, bl.ss = sei_1a, blp.ss = sei_1b, blf.ss = sei_1c)
fu.ss <- select(fu, resp_id, fu.ss = sei_1a, fup.ss = sei_1b, fuf.ss = sei_1c)
ss <- full_join(bl.ss, fu.ss, by = "resp_id")
full <- left_join(full, ss, by = "resp_id")

#####Identify as a Leader#####
bl.lead <- select(bl, resp_id, bl.lead = sei_13) 
pw.lead <- select(pw, resp_id, pw.lead = sei_13) 
fu.lead <- select(fu, resp_id, fu.lead = sei_13) 
lead <- full_join(bl.lead, fu.lead, by = "resp_id")
lead <- left_join(lead, pw.lead, by = "resp_id")
full <- left_join(full, lead, by = "resp_id")

#####Connectedness to Workshop Group#####
bl.tg <- select(bl, resp_id, bl.tg = sei_2) 
pw.tg <- select(pw, resp_id, pw.tg = sei_2) 
fu.tg <- select(fu, resp_id, fu.tg = sei_2) 
tg <- full_join(bl.tg, fu.tg, by = "resp_id")
tg <- left_join(tg, pw.tg, by = "resp_id")
full <- left_join(full, tg, by = "resp_id")

#####Aspirations for the Future#####
bl.a <- select(bl, resp_id, bl.rm = sei_14, bl.rmasp = sei_14a, sei_14b, bl.goal = sei_18a)
pw.a <- select(pw, resp_id, pw.rm = sei_14, pw.rmasp = sei_14a, sei_14b, pw.goal = sei_18ap)
fu.a <- select(fu, resp_id, fu.rm = sei_14, fu.rmasp = sei_14a, sei_14b, fu.goal = sei_18a, com.goal)

bl.a <- mutate(bl.a, bl.rmtime = if_else(sei_14b=="More than 10 years", 10, 
                                         if_else(sei_14b=="5-10 years", 7.5,
                                                 if_else(sei_14b=="3-5 years", 4,
                                                         if_else(sei_14b=="1-3 years",2,
                                                                 if_else(sei_14b=="6-12 months",0.75,
                                                                         if_else(sei_14b=="Less than 6 months",0.25,as.double(NA))))))),
               bl.imp = if_else(sei_14b=="It is impossible",1,0))

pw.a <- mutate(pw.a, pw.rmtime = if_else(sei_14b=="More than 10 years", 10, 
                                         if_else(sei_14b=="5-10 years", 7.5,
                                                 if_else(sei_14b=="3-5 years", 4,
                                                         if_else(sei_14b=="1-3 years",2,
                                                                 if_else(sei_14b=="6-12 months",0.75,
                                                                         if_else(sei_14b=="Less than 6 months",0.25,as.double(NA))))))),
               pw.imp = if_else(sei_14b=="It is impossible",1,0))

fu.a <- mutate(fu.a, fu.rmtime = if_else(sei_14b=="More than 10 years", 10, 
                                         if_else(sei_14b=="5-10 years", 7.5,
                                                 if_else(sei_14b=="3-5 years", 4,
                                                         if_else(sei_14b=="1-3 years",2,
                                                                 if_else(sei_14b=="6-12 months",0.75,
                                                                         if_else(sei_14b=="Less than 6 months",0.25,as.double(NA))))))),
               fu.imp = if_else(sei_14b=="It is impossible",1,0))

bl.a <- select(bl.a, resp_id, bl.rm, bl.rmasp, bl.rmtime, bl.imp, bl.goal)
pw.a <- select(pw.a, resp_id, pw.rm, pw.rmasp, pw.rmtime, pw.imp, pw.goal)
fu.a <- select(fu.a, resp_id, fu.rm, fu.rmasp, fu.rmtime, fu.imp, fu.goal, com.goal)
asp <- full_join(bl.a, fu.a, by = "resp_id")
asp <- left_join(asp, pw.a, by = "resp_id")

full <- left_join(full, asp, by = "resp_id")

#####Marginal Utility of Expenditure####
#Create datasets to calculate MUE at baseline and follow-up
bl.mue <- select(bl, resp_id, men = d_27a, women = d_27b, girls = d_27c, boys = d_27d, hhsize = d_26,
                 rice, onion, avocado, milk, soya, beans, fish, oil,
                 snacks, swpot, cabbage, cassava, sweets, mflour, tomato, beef, eggplant, banana)
bl.mue[] <- lapply(bl.mue, function(x) as.numeric(as.character(x)))

fu.mue <- select(fu, resp_id, men, women, girls, boys, hhsize, rice, onion, avocado, milk, 
                 soya, beans, fish, oil, snacks, swpot, cabbage, cassava, sweets, mflour, tomato, 
                 beef, eggplant, banana)
fu.mue[] <- lapply(fu.mue, function(x) as.numeric(as.character(x)))

bl.mue <- mutate(bl.mue, j = resp_id, t = 0, m = 0, lhhsize = log(hhsize))
fu.mue <- mutate(fu.mue, j = resp_id, t = 1, m = 0, lhhsize = log(hhsize))

#Create z matrix of household demographic characteristics
bl.z <- select(bl.mue, j, t, m, lhhsize, men, women, boys, girls)
fu.z <- select(fu.mue, j, t, m, lhhsize, men, women, boys, girls)
z <- rbind(bl.z, fu.z)
z <- filter(z, !(j==335 & t==1 & lhhsize < 2) & !(j==39 & t==1 & lhhsize < 1.9))

write.csv(z, "z_matrix.csv")

#Create y matrix of log expenditures
bl.mue <- mutate(bl.mue, lrice = log(rice), lonion = log(onion), lavocado = log(avocado), lmilk = log(milk),
                 lsoya = log(soya), lbeans = log(beans), lfish = log(fish), loil = log(oil),
                 lsnacks = log(snacks), lswpot = log(swpot), lcabbage = log(cabbage), lcassava = log(cassava),
                 lsweets = log(sweets), lmflour = log(mflour), ltomato = log(tomato), lbeef = log(beef),
                 leggplant = log(eggplant), lbanana = log(banana)) %>%
  select(j, t, m, lrice, lonion, lavocado, lmilk, lsoya, lbeans, lfish, loil, lsnacks, lswpot, lcabbage,
         lcassava, lsweets, lmflour, ltomato, lbeef, leggplant, lbanana)

fu.mue <- mutate(fu.mue, lrice = log(rice), lonion = log(onion), lavocado = log(avocado), lmilk = log(milk),
                 lsoya = log(soya), lbeans = log(beans), lfish = log(fish), loil = log(oil),
                 lsnacks = log(snacks), lswpot = log(swpot), lcabbage = log(cabbage), lcassava = log(cassava),
                 lsweets = log(sweets), lmflour = log(mflour), ltomato = log(tomato), lbeef = log(beef),
                 leggplant = log(eggplant), lbanana = log(banana)) %>%
  select(j, t, m, lrice, lonion, lavocado, lmilk, lsoya, lbeans, lfish, loil, lsnacks, lswpot, lcabbage,
         lcassava, lsweets, lmflour, ltomato, lbeef, leggplant, lbanana)

y <- rbind(bl.mue, fu.mue)
y[y==-Inf] <- NA #remove all negative infinities (0 consumption that was logged)
y <- filter(y, !(j==335 & t==1 & lrice < 7) & !(j==39 & t==1 & lrice<8))
write.csv(y, "y_matrix.csv")

#Read in MUE (calculated in python using a google colab notebook: https://colab.research.google.com/drive/1Hp-S9cUpCmzJQq3DjNv7yhiTI6KMKIdt?usp=sharing)
mue <- read.csv("IMUE.csv")

bl.mue <- filter(mue, t==0) %>% select(resp_id = j, bl.mue = loglambdas)
fu.mue <- filter(mue, t==1) %>% select(resp_id = j, fu.mue = loglambdas)

mue <- left_join(bl.mue, fu.mue, by = "resp_id")
full <- left_join(full, mue, by = "resp_id")

#####Consumption expenditures#####
bl.exp <- select(bl, resp_id, men = d_27a, women = d_27b, girls = d_27c, boys = d_27d, hhsize = d_26,
                 rice, onion, avocado, milk, soya, beans, fish, oil,
                 snacks, swpot, cabbage, cassava, sweets, mflour, tomato, beef, eggplant, banana)
bl.exp[] <- lapply(bl.exp, function(x) as.numeric(as.character(x)))

#Clean a few observations with unreasonably high responses that didn't get caught during data checks
bl.exp$rice[which(bl.exp$rice==30000)] <- 3000
bl.exp$swpot[which(bl.exp$swpot==14000)] <- 1400
bl.exp$swpot[which(bl.exp$swpot==15000)] <- 1500
bl.exp$beef[which(bl.exp$beef==37500)] <- 3750

bl.exp[] <- lapply(bl.exp, function(x) as.numeric(as.character(x)))
bl.exp <- mutate(bl.exp, total.exp = rice + onion + avocado + milk + soya + beans + fish + oil +
                   snacks + swpot + cabbage + cassava + sweets + mflour + tomato + beef + eggplant + banana,
                 exp.pc = total.exp/hhsize,
                 log.total.exp = log(total.exp),
                 log.exp.pc = log(exp.pc)) %>% select(resp_id, bl.total.exp = total.exp, bl.exp.pc = exp.pc, 
                                                      bl.log.total.exp = log.total.exp, bl.log.exp.pc = log.exp.pc)

fu.exp <- select(fu, resp_id, men, women, girls, boys, hhsize, rice, onion, avocado, milk, 
                 soya, beans, fish, oil, snacks, swpot, cabbage, cassava, sweets, mflour, tomato, 
                 beef, eggplant, banana)
fu.exp[] <- lapply(fu.exp, function(x) as.numeric(as.character(x)))

#Clean a few observations with unreasonably high responses that didn't get caught during data checks
fu.exp$rice[which(fu.exp$rice==20000)] <- 2000
fu.exp$rice[which(fu.exp$rice==18000)] <- 1800
fu.exp$rice[which(fu.exp$rice==16000)] <- 1600
fu.exp$beans[which(fu.exp$beans==15000)] <- 1500
fu.exp$swpot[which(fu.exp$swpot==16000)] <- 1600
fu.exp$mflour[which(fu.exp$mflour==12600)] <- 1260
fu.exp$mflour[which(fu.exp$mflour==12000)] <- 1200

fu.exp <- mutate(fu.exp, total.exp = rice + onion + avocado + milk + soya + beans + fish + oil +
                   snacks + swpot + cabbage + cassava + sweets + mflour + tomato + beef + eggplant + banana,
                 exp.pc = total.exp/hhsize,
                 log.total.exp = log(total.exp),
                 log.exp.pc = log(exp.pc))  %>% select(resp_id, fu.total.exp = total.exp, fu.exp.pc = exp.pc, 
                                                       fu.log.total.exp = log.total.exp, fu.log.exp.pc = log.exp.pc)

exp <- left_join(bl.exp, fu.exp, by = "resp_id")
full <- left_join(full, exp, by = "resp_id")

#####Income and Work#####
bl.minc <- select(bl, resp_id, e_21, e_22, e_22_int, e_19, e_19_other,
                  e_20, e_20_other) %>%
  mutate(bl.income = if_else(e_22_int=="RWF 0 - 4999", 2500,
                             if_else(e_22_int=="RWF 10,000 - 14,999", 12500,
                                     if_else(e_22_int=="RWF 20,000 - 24,999", 22500,
                                             if_else(e_22_int=="RWF 75,000 - 79,999",77500,
                                                     if_else(e_22_int=="Above RWF 80,000",80000,
                                                             if_else(e_22_int=="Don't know",as.double(NA),
                                                                     if_else(e_21=="No",0,as.double(e_22)))))))),
         bl.work_lastmonth = e_21,
         bl.ownfarm = if_else(e_19=="Farming – own/family farm",1,0),
         bl.paidfarm = if_else(e_19=="Farming – paid worker on other’s farm",1,0),
         bl.selfwork = if_else(e_19=="Shopkeeping – own shop/family shop"|e_19=="Other [specify]:___________________",1,0),
         bl.selfemp = if_else(e_20=="Self-employed (work for yourself)",1,0),
         bl.paidemp = if_else(e_20=="In paid employment (work for someone else)",1,0),
         bl.cantwork = if_else(e_20=="Unable to work due to sickness or ill-health",1,0),
         bl.searchwork = if_else(e_20=="Looking for work",1,0),
         bl.home = if_else(e_20=="Looking after the home",1,0),
         bl.notinlf = if_else(e_20=="Not working and not looking for work",1,0),
         bl.ret = if_else(e_20=="Retired",1,0))

fu.minc <- select(fu, resp_id, e_21, e_22, e_22_int, e_19, e_19_other, e_20, e_20_other)
fu.minc <- mutate(fu.minc, fu.income = if_else(e_21=="No",0,as.double(e_22)),
                  fu.work_lastmonth = e_21,
                  fu.ownfarm = if_else(e_19=="Farming – own/family farm",1,0),
                  fu.paidfarm = if_else(e_19=="Farming – paid worker on other’s farm",1,0),
                  fu.selfwork = if_else(e_19=="Shopkeeping – own shop/family shop"|e_19=="Other",1,0),
                  fu.selfemp = if_else(e_20=="Self-employed (work for yourself)",1,0),
                  fu.paidemp = if_else(e_20=="In paid employment (work for someone else)",1,0),
                  fu.cantwork = if_else(e_20=="Unable to work due to sickness or ill-health"|e_20=="Other",1,0),
                  fu.searchwork = if_else(e_20=="Looking for work",1,0),
                  fu.home = if_else(e_20=="Looking after the home",1,0),
                  fu.notinlf = if_else(e_20=="Student",1,0))

bl.minc <- select(bl.minc, -e_21, -e_22, -e_22_int, -e_19, -e_19_other, -e_20, -e_20_other)
bl.minc$bl.income[which(is.na(bl.minc$bl.income))] <- 0
fu.minc <- select(fu.minc, -e_21, -e_22, -e_19, -e_19_other, -e_20, -e_20_other)
fu.minc$fu.income[which(is.na(fu.minc$fu.income))] <- 0
minc <- full_join(bl.minc, fu.minc, by = "resp_id")

minc[, colnames(minc)] <- lapply(colnames(minc), function(x) as.numeric(minc[[x]]))

full <- left_join(full, minc, by = "resp_id")

full <- mutate(full, ihs.fu = log(fu.income + sqrt(fu.income^2 + 1)),
               ihs.bl = log(bl.income + sqrt(bl.income^2 + 1)))

#####Speaking up####
bl.su <- select(bl, resp_id, bl.su = sei_8a)
fu.su <- select(fu, resp_id, fu.su = sei_8a)
su <- full_join(bl.su, fu.su, by = "resp_id")
full <- left_join(full, su, by = "resp_id")
#####Remove extraneous dataframes and create treatment variables#####
#rm(list=setdiff(ls(), "full","bl","pw","fu"))

full <- mutate(full, SFL = if_else(workshop_name=="Storytelling for Leadership",1,0),
               PD = if_else(workshop_name=="Professional Development",1,0),
               tmt = if_else(workshop_name!="Control",1,0))

#####Merge in post workshop data#####
pw.merge <- pw[,c(7,40:57)]
full <- left_join(full, pw.merge, by = "resp_id")

#####Construct Anderson (2008) index variables of all outcomes in each pre-specified family####
full <- mutate(full, bl.mue2 = bl.mue*(-1), fu.mue2 = fu.mue*(-1))
pps.bl.means <- filter(full, workshop_name=="Control") %>%
  select(bl.saew, bl.conhvew, bl.conew, bl.comew, bl.psew, bl.lead.ind, bl.nlead,
         bl.ss, blp.ss, blf.ss, bl.lead, bl.tg, bl.rm, bl.rmasp, bl.rmtime, bl.rew, bl.vew, bl.goal,
         bl.mue2, bl.log.exp.pc, bl.work_lastmonth, bl.income) %>% 
  summarise_all(list(mean), na.rm = T)
pps.bl.sds <- filter(full, workshop_name=="Control") %>%
  select(bl.saew, bl.conhvew, bl.conew, bl.comew, bl.psew, bl.lead.ind, bl.nlead,
         bl.ss, blp.ss, blf.ss, bl.lead, bl.tg, bl.rm, bl.rmasp, bl.rmtime, bl.rew, bl.vew, bl.goal,
         bl.mue2, bl.log.exp.pc, bl.work_lastmonth, bl.income) %>% 
  summarise_all(list(sd), na.rm = T)

pps.fu.means <- filter(full, workshop_name=="Control") %>%
  select(fu.saew, fu.conhvew, fu.conew, fu.comew, fu.psew, fu.lead.ind, fu.nlead,
         fu.ss, fup.ss, fuf.ss, fu.lead, fu.tg, fu.rm, fu.rmasp, fu.rmtime, fu.rew, fu.vew, fu.goal, com.goal,
         fu.mue2, fu.log.exp.pc, fu.work_lastmonth, fu.income) %>% 
  summarise_all(list(mean), na.rm=T)
pps.fu.sds <- filter(full, workshop_name=="Control") %>%
  select(fu.saew, fu.conhvew, fu.conew, fu.comew, fu.psew, fu.lead.ind, fu.nlead,
         fu.ss, fup.ss, fuf.ss, fu.lead, fu.tg, fu.rm, fu.rmasp, fu.rmtime, fu.rew, fu.vew, fu.goal, com.goal,
         fu.mue2, fu.log.exp.pc, fu.work_lastmonth, fu.income) %>% 
  summarise_all(list(sd), na.rm=T)

# Demean and convert to SD units
full <- mutate(full, blsa2 = (bl.saew - pps.bl.means[1,1])/pps.bl.sds[1,1],
               blconhv2 = (bl.conhvew - pps.bl.means[1,2])/pps.bl.sds[1,2],
               blcon2 = (bl.conew - pps.bl.means[1,3])/pps.bl.sds[1,3],
               blcom2 = (bl.comew - pps.bl.means[1,4])/pps.bl.sds[1,4],
               blps2 = (bl.psew - pps.bl.means[1,5])/pps.bl.sds[1,5],
               blleadind2 = (bl.lead.ind - pps.bl.means[1,6])/pps.bl.sds[1,6],
               blnlead2 = (bl.nlead - pps.bl.means[1,7])/pps.bl.sds[1,7],
               blss2 = (bl.ss - pps.bl.means[1,8])/pps.bl.sds[1,8],
               blpss2 = (blp.ss - pps.bl.means[1,9])/pps.bl.sds[1,9],
               blfss2 = (blf.ss - pps.bl.means[1,10])/pps.bl.sds[1,10],
               bllead2 = (bl.lead - pps.bl.means[1,11])/pps.bl.sds[1,11],
               bltg2 = (bl.tg - pps.bl.means[1,12])/pps.bl.sds[1,12],
               blrm2 = (bl.rm - pps.bl.means[1,13])/pps.bl.sds[1,13],
               blrmasp2 = (bl.rmasp - pps.bl.means[1,14])/pps.bl.sds[1,14],
               blrmtime2 = -1*(bl.rmtime - pps.bl.means[1,15])/pps.bl.sds[1,15],
               blr2 = (bl.rew - pps.bl.means[1,16])/pps.bl.sds[1,16],
               blv2 = (bl.vew - pps.bl.means[1,17])/pps.bl.sds[1,17],
               blgoal2 = (bl.goal - pps.bl.means[1,18])/pps.bl.sds[1,18],
               blmue2 = (bl.mue2 - pps.bl.means[1,19])/pps.bl.sds[1,19],
               blexp = (bl.log.exp.pc - pps.bl.means[1,20])/pps.bl.sds[1,20],
               blwork2 = (bl.work_lastmonth - pps.bl.means[1,21])/pps.bl.sds[1,21],
               blinc2 = (bl.income - pps.bl.means[1,22])/pps.bl.sds[1,22],
               fusa2 = (fu.saew - pps.fu.means[1,1])/pps.fu.sds[1,1],
               fuconhv2 = (fu.conhvew - pps.fu.means[1,2])/pps.fu.sds[1,2],
               fucon2 = (fu.conew - pps.fu.means[1,3])/pps.fu.sds[1,3],
               fucom2 = (fu.comew - pps.fu.means[1,4])/pps.fu.sds[1,4],
               fups2 = (fu.psew - pps.fu.means[1,5])/pps.fu.sds[1,5],
               fuleadind2 = (fu.lead.ind - pps.fu.means[1,6])/pps.fu.sds[1,6],
               funlead2 = (fu.nlead - pps.fu.means[1,7])/pps.fu.sds[1,7],
               fuss2 = (fu.ss - pps.fu.means[1,8])/pps.fu.sds[1,8],
               fupss2 = (fup.ss - pps.fu.means[1,9])/pps.fu.sds[1,9],
               fufss2 = (fuf.ss - pps.fu.means[1,10])/pps.fu.sds[1,10],
               fulead2 = (fu.lead - pps.fu.means[1,11])/pps.fu.sds[1,11],
               futg2 = (fu.tg - pps.fu.means[1,12])/pps.fu.sds[1,12],
               furm2 = (fu.rm - pps.fu.means[1,13])/pps.fu.sds[1,13],
               furmasp2 = (fu.rmasp - pps.fu.means[1,14])/pps.fu.sds[1,14],
               furmtime2 = -1*(fu.rmtime - pps.fu.means[1,15])/pps.fu.sds[1,15],
               fur2 = (fu.rew - pps.fu.means[1,16])/pps.fu.sds[1,16],
               fuv2 = (fu.vew - pps.fu.means[1,17])/pps.fu.sds[1,17],
               fugoal2 = (fu.goal - pps.fu.means[1,18])/pps.fu.sds[1,18],
               fucomgoal2 = (com.goal - pps.fu.means[1,19])/pps.fu.sds[1,19],
               fumue2 = (fu.mue2 - pps.fu.means[1,20])/pps.fu.sds[1,20],
               fuexp = (fu.log.exp.pc - pps.fu.means[1,21])/pps.fu.sds[1,21],
               fuwork2 = (fu.work_lastmonth - pps.fu.means[1,22])/pps.fu.sds[1,22],
               fuinc2 = (fu.income - pps.fu.means[1,23])/pps.fu.sds[1,23])

#Calculate the variance/covariance matrix for each survey round and each family of outcomes
bl.pps.vars <- select(full, blsa2, blconhv2, blcon2, blcom2, blps2, blleadind2, blnlead2) 
bl.pps.vars <- bl.pps.vars[complete.cases(bl.pps.vars),]
fu.pps.vars <- select(full, fusa2, fuconhv2, fucon2, fucom2, fups2, fuleadind2, funlead2)
fu.pps.vars <- fu.pps.vars[complete.cases(fu.pps.vars),]

#Note that I'm going to leave out aspires to be role model and time to be role model because there are a relatively large number of missing values
bl.sps.vars <- select(full, blss2, blpss2, blfss2, bllead2, bltg2, blrm2, blr2, blv2, blgoal2) 
bl.sps.vars <- bl.sps.vars[complete.cases(bl.sps.vars),]
fu.sps.vars <- select(full, fuss2, fupss2, fufss2, fulead2, futg2, furm2, fur2, fuv2, fugoal2, fucomgoal2)
fu.sps.vars <- fu.sps.vars[complete.cases(fu.sps.vars),]

bl.e.vars <- select(full, blmue2, blexp, blwork2, blinc2) #%>% mutate(blmue = blmue2*-1) %>% select(-blmue2)
bl.e.vars <- bl.e.vars[complete.cases(bl.e.vars),]
fu.e.vars <- select(full, fumue2, fuexp, fuwork2, fuinc2)  #%>% mutate(fumue = fumue2*-1) %>% select(-fumue2)
fu.e.vars <- fu.e.vars[complete.cases(fu.e.vars),]

bl.pps.vcov <- cov(bl.pps.vars)
fu.pps.vcov <- cov(fu.pps.vars)

bl.sps.vcov <- cov(bl.sps.vars)
fu.sps.vcov <- cov(fu.sps.vars)

bl.e.vcov <- cov(bl.e.vars)
fu.e.vcov <- cov(fu.e.vars)

pps.ones <- c(1,1,1,1,1,1,1)
bl.pps.coef <- solve(t(pps.ones)%*%solve(bl.pps.vcov)%*%pps.ones)
fu.pps.coef <- solve(t(pps.ones)%*%solve(fu.pps.vcov)%*%pps.ones)

sps.ones.b <- c(1,1,1,1,1,1,1,1,1)
sps.ones.f <- c(1,1,1,1,1,1,1,1,1,1)
bl.sps.coef <- solve(t(sps.ones.b)%*%solve(bl.sps.vcov)%*%sps.ones.b)
fu.sps.coef <- solve(t(sps.ones.f)%*%solve(fu.sps.vcov)%*%sps.ones.f)

e.ones <- c(1,1,1,1)
bl.e.coef <- solve(t(e.ones)%*%solve(bl.e.vcov)%*%e.ones)
fu.e.coef <- solve(t(e.ones)%*%solve(fu.e.vcov)%*%e.ones)

full$bl.pps.ind <- NA
full$fu.pps.ind <- NA
full$bl.sps.ind <- NA
full$fu.sps.ind <- NA
full$bl.e.ind <- NA
full$fu.e.ind <- NA

for (i in 1:length(full$resp_id)){
  bl.pps.outcomes <- as.matrix(full[i,c(120:126)])
  fu.pps.outcomes <- as.matrix(full[i,c(142:148)])
  bl.sps.outcomes <- as.matrix(full[i,c(127:132,135:137)])
  fu.sps.outcomes <- as.matrix(full[i,c(149:154,157:160)])
  bl.e.outcomes <- as.matrix(full[i,c(138:141)])
  fu.e.outcomes <- as.matrix(full[i,c(161:164)])
  
  full$bl.pps.ind[i] <- bl.pps.coef%*%t(pps.ones)%*%solve(bl.pps.vcov)%*%t(bl.pps.outcomes)
  full$fu.pps.ind[i] <- fu.pps.coef%*%t(pps.ones)%*%solve(fu.pps.vcov)%*%t(fu.pps.outcomes)
  
  full$bl.sps.ind[i] <- bl.sps.coef%*%t(sps.ones.b)%*%solve(bl.sps.vcov)%*%t(bl.sps.outcomes)
  full$fu.sps.ind[i] <- fu.sps.coef%*%t(sps.ones.f)%*%solve(fu.sps.vcov)%*%t(fu.sps.outcomes)
  
  full$bl.e.ind[i] <- bl.e.coef%*%t(e.ones)%*%solve(bl.e.vcov)%*%t(bl.e.outcomes)
  full$fu.e.ind[i] <- fu.e.coef%*%t(e.ones)%*%solve(fu.e.vcov)%*%t(fu.e.outcomes)
}

#####Merge in data on workshops and cooperatives####
full <- left_join(full, coop, by = c("resp_id" = "ID"))
ws <- select(wksp, -Coop)
full <- left_join(full, ws, by = "resp_id")
##########################################################################################
#Pre-Specified Regressions - Primary Socio-Emotional Outcomes
##########################################################################################
##NOTE: for each outcome where we have baseline observations, first estimate the ANCOVA 
##then estimate omitting the baseline control for robustness

#Self-advocacy
sa.reg.ew <- felm(fu.saew ~ SFL + PD + bl.saew|0|0, data = full)
summary(sa.reg.ew, robust = TRUE)
linearHypothesis(sa.reg.ew, "SFL  = PD")

sa.reg.ew2 <- felm(fu.saew ~ SFL + PD|0|0, data = full)
summary(sa.reg.ew2, robust = TRUE)
linearHypothesis(sa.reg.ew2, "SFL  = PD")

#Problem Solving
ps.reg.ew <- felm(fu.psew ~ SFL + PD + bl.psew|0|0, data = full)
summary(ps.reg.ew, robust = TRUE)
linearHypothesis(ps.reg.ew, "SFL  = PD")

ps.reg.ew2 <- felm(fu.psew ~ SFL + PD|0|0, data = full)
summary(ps.reg.ew2, robust = TRUE)
linearHypothesis(ps.reg.ew2, "SFL  = PD")

#No. Leadership positions
nlead.reg <- felm(fu.nlead ~ SFL + PD + bl.nlead|0|0, data = full)
summary(nlead.reg, robust = TRUE)
linearHypothesis(nlead.reg, "SFL  = PD")

nlead.reg2 <- felm(fu.nlead ~ SFL + PD|0|0, data = full)
summary(nlead.reg2, robust = TRUE)
linearHypothesis(nlead.reg2, "SFL  = PD")

#Connections to village/home
conhv.reg.ew <- felm(fu.conhvew ~ SFL + PD + bl.conhvew|0|0, data = full)
summary(conhv.reg.ew, robust = TRUE)
linearHypothesis(conhv.reg.ew, "SFL  = PD")

conhv.reg.ew2 <- felm(fu.conhvew ~ SFL + PD|0|0, data = full)
summary(conhv.reg.ew2, robust = TRUE)
linearHypothesis(conhv.reg.ew2, "SFL  = PD")

#Connections to cooperatives 
con.reg.ew <- felm(fu.conew ~ SFL + PD + bl.conew|0|0, data = full)
summary(con.reg.ew, robust = TRUE)
linearHypothesis(con.reg.ew, "SFL  = PD")

con.reg.ew2 <- felm(fu.conew ~ SFL + PD|0|0, data = full)
summary(con.reg.ew2, robust = TRUE)
linearHypothesis(con.reg.ew2, "SFL  = PD")

#Communication
com.reg.ew <- felm(fu.comew ~ SFL + PD + bl.comew|0|0, data = full)
summary(com.reg.ew, robust = TRUE)
linearHypothesis(com.reg.ew, "SFL  = PD")

com.reg.ew2 <- felm(fu.comew ~ SFL + PD|0|0, data = full)
summary(com.reg.ew2, robust = TRUE)
linearHypothesis(com.reg.ew2, "SFL  = PD")

#Any leadership position
lead.ind.reg <- felm(fu.lead.ind ~ SFL + PD + bl.lead.ind|0|0, data = full)
summary(lead.ind.reg, robust = TRUE)
linearHypothesis(lead.ind.reg, "SFL  = PD")

lead.ind.reg2 <- felm(fu.lead.ind ~ SFL + PD |0|0, data = full)
summary(lead.ind.reg2, robust = TRUE)
linearHypothesis(lead.ind.reg2, "SFL  = PD")

##########################################################################################
#Pre-Specified Regressions - Secondary Socio-Emotional Outcomes
##########################################################################################
#Social status - past, future, current
ssp.reg <- felm(fup.ss ~ SFL + PD + blp.ss|0|0, data = full)
summary(ssp.reg, robust = TRUE)
linearHypothesis(ssp.reg, "SFL  = PD")

ssp.reg.2 <- felm(fup.ss ~ SFL + PD |0|0, data = full)
summary(ssp.reg.2, robust = TRUE)
linearHypothesis(ssp.reg.2, "SFL  = PD")

ssf.reg <- felm(fuf.ss ~ SFL + PD + blf.ss|0|0, data = full)
summary(ssf.reg, robust = TRUE)
linearHypothesis(ssf.reg, "SFL  = PD")

ssf.reg.2 <- felm(fuf.ss ~ SFL + PD|0|0, data = full)
summary(ssf.reg.2, robust = TRUE)
linearHypothesis(ssf.reg.2, "SFL  = PD")

ss.reg <- felm(fu.ss ~ SFL + PD + bl.ss|0|0, data = full)
summary(ss.reg, robust = TRUE)
linearHypothesis(ss.reg, "SFL  = PD")

ss.reg.2 <- felm(fu.ss ~ SFL + PD |0|0, data = full)
summary(ss.reg.2, robust = TRUE)
linearHypothesis(ss.reg.2, "SFL  = PD")

#Self-value
val.reg.ew <- felm(fu.vew ~ SFL + PD + bl.vew|0|0, data = full)
summary(val.reg.ew, robust = TRUE)
linearHypothesis(val.reg.ew, "SFL  = PD")

val.reg.ew2 <- felm(fu.vew ~ SFL + PD |0|0, data = full)
summary(val.reg.ew2, robust = TRUE)
linearHypothesis(val.reg.ew2, "SFL  = PD")

#Resilience
res.reg.ew <- felm(fu.rew ~ SFL + PD + bl.rew|0|0, data = full)
summary(res.reg.ew, robust = TRUE)
linearHypothesis(res.reg.ew, "SFL  = PD")

res.reg.ew2 <- felm(fu.rew ~ SFL + PD|0|0, data = full)
summary(res.reg.ew2, robust = TRUE)
linearHypothesis(res.reg.ew2, "SFL  = PD")

#Identify as a leader
lead.reg <- felm(fu.lead ~ SFL + PD + bl.lead|0|0, data = full)
summary(lead.reg, robust = TRUE)
linearHypothesis(lead.reg, "SFL  = PD")

lead.reg2 <- felm(fu.lead ~ SFL + PD |0|0, data = full)
summary(lead.reg2, robust = TRUE)
linearHypothesis(lead.reg2, "SFL  = PD")

#Connections to training group
tg.reg <- felm(fu.tg ~ SFL + PD + bl.tg|0|0, data = full)
summary(tg.reg, robust = TRUE)
linearHypothesis(tg.reg, "SFL  = PD")

tg.reg2 <- felm(fu.tg ~ SFL + PD|0|0, data = full)
summary(tg.reg2, robust = TRUE)
linearHypothesis(tg.reg2, "SFL  = PD")

#Time to be role model
rmtime.reg <- felm(fu.rmtime ~ SFL + PD + bl.rmtime|0|0, data = full)
summary(rmtime.reg, robust = TRUE)
linearHypothesis(rmtime.reg, "SFL  = PD")

rmtime.reg.2 <- felm(fu.rmtime ~ SFL + PD |0|0, data = full)
summary(rmtime.reg.2, robust = TRUE)
linearHypothesis(rmtime.reg.2, "SFL  = PD")

#Has role model
rm.reg <- felm(fu.rm ~ SFL + PD + bl.rm|0|0, data = full)
summary(rm.reg, robust = TRUE)
linearHypothesis(rm.reg, "SFL  = PD")

rm.reg.2 <- felm(fu.rm ~ SFL + PD |0|0, data = full)
summary(rm.reg.2, robust = TRUE)
linearHypothesis(rm.reg.2, "SFL  = PD")

#Has goal
goal.reg <- felm(fu.goal ~ SFL + PD + bl.goal|0|0, data = full)
summary(goal.reg, robust = TRUE)
linearHypothesis(goal.reg, "SFL  = PD")

goal.reg.2 <- felm(fu.goal ~ SFL + PD |0|0, data = full)
summary(goal.reg.2, robust = TRUE)
linearHypothesis(goal.reg.2, "SFL  = PD")

#Aspires to be role model
rmasp.reg <- felm(fu.rmasp ~ SFL + PD + bl.rmasp|0|0, data = full)
summary(rmasp.reg, robust = TRUE)
linearHypothesis(rmasp.reg, "SFL  = PD")

rmasp.reg.2 <- felm(fu.rmasp ~ SFL + PD |0|0, data = full)
summary(rmasp.reg.2, robust = TRUE)
linearHypothesis(rmasp.reg.2, "SFL  = PD")

#Achieved last goal
com.goal.reg <- felm(com.goal ~ SFL + PD|0|0, data = full)
summary(com.goal.reg, robust = TRUE)
linearHypothesis(com.goal.reg, "SFL  = PD")

##########################################################################################
#Pre-Specified Regressions - Economic Outcomes
##########################################################################################
#Income
inc <- felm(fu.income ~ SFL + PD + bl.income |0|0, data = full)
summary(inc, robust = TRUE)
linearHypothesis(inc, "SFL = PD")

inc2 <- felm(fu.income ~ SFL + PD |0|0, data = full)
summary(inc2, robust = TRUE)
linearHypothesis(inc, "SFL = PD")

#Earn any income
work.reg <- felm(fu.work_lastmonth ~ SFL + PD + bl.work_lastmonth|0|0, data = full)
summary(work.reg, robust = TRUE)
linearHypothesis(work.reg, "SFL = PD")

work.reg2 <- felm(fu.work_lastmonth ~ SFL + PD |0|0, data = full)
summary(work.reg2, robust = TRUE)
linearHypothesis(work.reg2, "SFL = PD")

#MUE
mue.reg <- felm(fu.mue ~ SFL + PD + bl.mue|0|0, data = full)
summary(mue.reg, robust = TRUE)
linearHypothesis(mue.reg, "SFL = PD")

mue.reg2 <- felm(fu.mue ~ SFL + PD |0|0, data = full)
summary(mue.reg, robust = TRUE)
linearHypothesis(mue.reg2, "SFL = PD")

#Not pre-specified, but requested by reviewers: consumption expenditures
exp.reg <- felm(fu.log.exp.pc ~ SFL + PD + bl.log.exp.pc|0|0, data = full)
summary(exp.reg, robust = TRUE)
linearHypothesis(exp.reg, "SFL = PD")

exp.reg2 <- felm(fu.log.total.exp ~ SFL + PD|0|0, data = full)
summary(exp.reg2, robust = TRUE)
linearHypothesis(exp.reg2, "SFL = PD")

##########################################################################################
#NOT Pre-Specified Regressions/responses to reviewer requests 
##########################################################################################
#####Probability of speaking up in a group situation####
su.reg <- felm(fu.su ~ SFL + PD + bl.su|0|0, data = full)
summary(su.reg, robust = TRUE)
linearHypothesis(su.reg, "SFL = PD")

#####Heterogeneity by age####
full <- mutate(full, over40 = if_else(age>=40,1,0),
               under40 = if_else(age<40,1,0))

#Primary psychosocial
sa.het <- felm(fu.saew ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.saew|0|0, data = full)
sa.het2 <- felm(fu.saew ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.saew|0|0, data = full)

ps.het <- felm(fu.psew ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.psew|0|0, data = full)
ps.het2 <- felm(fu.psew ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.psew|0|0, data = full)

nl.het <- felm(fu.nlead ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.nlead|0|0, data = full)
nl.het2 <- felm(fu.nlead ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.nlead|0|0, data = full)

li.het <- felm(fu.lead.ind ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.lead.ind|0|0, data = full)
li.het2 <- felm(fu.lead.ind ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.lead.ind|0|0, data = full)

chv.het <- felm(fu.conhvew ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.conhvew|0|0, data = full)
chv.het2 <- felm(fu.conhvew ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.conhvew|0|0, data = full)

cc.het <- felm(fu.conew ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.conew|0|0, data = full)
cc.het2 <- felm(fu.conew ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.conew|0|0, data = full)

com.het <- felm(fu.comew ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.comew|0|0, data = full)
com.het2 <- felm(fu.comew ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.comew|0|0, data = full)

pps.het <- felm(fu.pps.ind ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.pps.ind|0|0, data = full)
pps.het2 <- felm(fu.pps.ind ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.pps.ind|0|0, data = full)

#Secondary psychosocial outcomes
pss.het <- felm(fup.ss ~ SFL + PD + over40 + SFL:over40 + PD:over40 + blp.ss|0|0, data = full)
pss.het2 <- felm(fup.ss ~ SFL + PD + under40 + SFL:under40 + PD:under40 + blp.ss|0|0, data = full)

fss.het <- felm(fuf.ss ~ SFL + PD + over40 + SFL:over40 + PD:over40 + blf.ss|0|0, data = full)
fss.het2 <- felm(fuf.ss ~ SFL + PD + under40 + SFL:under40 + PD:under40 + blf.ss|0|0, data = full)

ss.het <- felm(fu.ss ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.ss|0|0, data = full)
ss.het2 <- felm(fu.ss ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.ss|0|0, data = full)

v.het <- felm(fu.vew ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.vew|0|0, data = full)
v.het2 <- felm(fu.vew ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.vew|0|0, data = full)

r.het <- felm(fu.rew ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.rew|0|0, data = full)
r.het2 <- felm(fu.rew ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.rew|0|0, data = full)

lid.het <- felm(fu.lead ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.lead|0|0, data = full)
lid.het2 <- felm(fu.lead ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.lead|0|0, data = full)

con.het <- felm(fu.tg~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.tg|0|0, data = full)
con.het2 <- felm(fu.tg ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.tg|0|0, data = full)

rmt.het <- felm(fu.rmtime ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.rmtime|0|0, data = full)
rmt.het2 <- felm(fu.rmtime ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.rmtime|0|0, data = full)

rm.het <- felm(fu.rm ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.rm|0|0, data = full)
rm.het2 <- felm(fu.rm ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.rm|0|0, data = full)

goal.het <- felm(fu.goal ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.goal|0|0, data = full)
goal.het2 <- felm(fu.goal ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.goal|0|0, data = full)

rma.het <- felm(fu.rmasp ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.rmasp|0|0, data = full)
rma.het2 <- felm(fu.rmasp ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.rmasp|0|0, data = full)

gc.het <- felm(com.goal ~ SFL + PD + over40 + SFL:over40 + PD:over40|0|0, data = full)
gc.het2 <- felm(com.goal ~ SFL + PD + under40 + SFL:under40 + PD:under40|0|0, data = full)

sps.het <- felm(fu.sps.ind ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.sps.ind|0|0, data = full)
sps.het2 <- felm(fu.sps.ind ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.sps.ind|0|0, data = full)

#Economic outcomes
mue.het <- felm(fu.mue ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.mue|0|0, data = full)
mue.het2 <- felm(fu.mue ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.mue|0|0, data = full)

i.het <- felm(fu.income ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.income|0|0, data = full)
i.het2 <- felm(fu.income ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.income|0|0, data = full)

c.het <- felm(fu.log.exp.pc ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.log.exp.pc|0|0, data = full)
c.het2 <- felm(fu.log.exp.pc ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.log.exp.pc|0|0, data = full)

ai.het <- felm(fu.work_lastmonth ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.work_lastmonth|0|0, data = full)
ai.het2 <- felm(fu.work_lastmonth ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.work_lastmonth|0|0, data = full)

e.het <- felm(fu.e.ind ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.e.ind|0|0, data = full)
e.het2 <- felm(fu.e.ind ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.e.ind|0|0, data = full)

#Speaking up
su.het <- felm(fu.su ~ SFL + PD + over40 + SFL:over40 + PD:over40 + bl.su|0|0, data = full)
su.het.2 <- felm(fu.su ~ SFL + PD + under40 + SFL:under40 + PD:under40 + bl.su|0|0, data = full)

#####Heterogeneity by positive cash earnings at baseline####
earn.pps <- felm(fu.pps.ind ~ SFL + PD + SFL:bl.work_lastmonth + 
                   PD:bl.work_lastmonth + bl.pps.ind + bl.work_lastmonth|0|0, data = full)
summary(earn.pps, robust = TRUE)
linearHypothesis(earn.pps, "SFL = PD")
linearHypothesis(earn.pps, "SFL:bl.work_lastmonth = PD:bl.work_lastmonth")

earn.sps <- felm(fu.sps.ind ~ SFL + PD + SFL:bl.work_lastmonth + 
                   PD:bl.work_lastmonth + bl.sps.ind + bl.work_lastmonth|0|0, data = full)
summary(earn.sps, robust = TRUE)
linearHypothesis(earn.sps, "SFL = PD")
linearHypothesis(earn.sps, "SFL:bl.work_lastmonth = PD:bl.work_lastmonth")

earn.e <- felm(fu.e.ind ~ SFL + PD + SFL:bl.work_lastmonth + 
                 PD:bl.work_lastmonth + bl.e.ind + bl.work_lastmonth|0|0, data = full)
summary(earn.e, robust = TRUE)
linearHypothesis(earn.e, "SFL = PD")
linearHypothesis(earn.e, "SFL:bl.work_lastmonth = PD:bl.work_lastmonth")

#####Heterogeneity by above median PPS####
median.pps <- median(full$bl.pps.ind)
full <- mutate(full, high.pps = if_else(bl.pps.ind>=median.pps,1,0),
               low.pps = if_else(bl.pps.ind<median.pps,1,0))

het.pps <- felm(fu.pps.ind ~ SFL + PD + SFL:high.pps + 
                  PD:high.pps + bl.pps.ind + high.pps|0|0, data = full)
het.pps2 <- felm(fu.pps.ind ~ SFL + PD + SFL:low.pps + 
                   PD:low.pps + bl.pps.ind + high.pps|0|0, data = full)
summary(het.pps, robust = TRUE)
linearHypothesis(het.pps, "SFL = PD")
linearHypothesis(het.pps, "SFL:high.pps = PD:high.pps")

het.sps <- felm(fu.pps.ind ~ SFL + PD + SFL:high.pps + 
                  PD:high.pps + bl.sps.ind + high.pps|0|0, data = full)
het.sps2 <- felm(fu.pps.ind ~ SFL + PD + SFL:low.pps + 
                   PD:low.pps + bl.sps.ind + high.pps|0|0, data = full)
summary(het.sps, robust = TRUE)
linearHypothesis(het.sps, "SFL = PD")
linearHypothesis(het.sps, "SFL:high.pps = PD:high.pps")

het.e <- felm(fu.e.ind ~ SFL + PD + SFL:high.pps + 
                PD:high.pps + bl.e.ind + high.pps|0|0, data = full)
het.e2 <- felm(fu.e.ind ~ SFL + PD + SFL:low.pps + 
                 PD:low.pps + bl.e.ind + high.pps|0|0, data = full)
summary(het.e, robust = TRUE)
linearHypothesis(het.e, "SFL = PD")
linearHypothesis(het.e, "SFL:high.pps = PD:high.pps")

names <- c(rep(c("SFL - Below Median", "PD - Below Median", 
                 "SFL - Above Median", "PD - Above Median"),3))
tmt <- c(rep(c("SFL","PD","SFL","PD"),3))
group <- c(rep(c("Below Median","Below Median","Above Median","Above Median"),3))
outcome <- c(rep("PPS Index",4),rep("SPS Index",4),rep("Econ Index",4))
coefs <- c(het.pps$coefficients[2], het.pps$coefficients[3],
           het.pps2$coefficients[2], het.pps2$coefficients[3],
           het.sps$coefficients[2], het.sps$coefficients[3],
           het.sps2$coefficients[2], het.sps2$coefficients[3],
           het.e$coefficients[2], het.e$coefficients[3],
           het.e2$coefficients[2], het.e2$coefficients[3])
SEs <- c(het.pps$rse[2], het.pps$rse[3],
         het.pps2$rse[2], het.pps2$rse[3],
         het.sps$rse[2], het.sps$rse[3],
         het.sps2$rse[2], het.sps2$rse[3],
         het.e$rse[2], het.e$rse[3],
         het.e2$rse[2], het.e2$rse[3])

het.pps.pdata <- as.data.frame(cbind(names, tmt, group, outcome, coefs, SEs))
het.pps.pdata$coefs <- as.numeric(het.pps.pdata$coefs)
het.pps.pdata$SEs <- as.numeric(het.pps.pdata$SEs)
het.pps.pdata <- mutate(het.pps.pdata, CI.High = coefs + 1.96*SEs,
                        CI.Low = coefs - 1.96*SEs)

#####Heterogeneity by above median SPS####
median.sps <- median(full$bl.sps.ind, na.rm=T)
full <- mutate(full, high.sps = if_else(bl.sps.ind>=median.sps,1,0),
               low.sps = if_else(bl.sps.ind<median.sps,1,0))

het2.pps <- felm(fu.pps.ind ~ SFL + PD + SFL:high.sps + 
                   PD:high.sps + bl.pps.ind + high.sps|0|0, data = full)
het2.pps2 <- felm(fu.pps.ind ~ SFL + PD + SFL:low.sps + 
                    PD:low.sps + bl.pps.ind + high.sps|0|0, data = full)
summary(het2.pps, robust = TRUE)
linearHypothesis(het2.pps, "SFL = PD")
linearHypothesis(het2.pps, "SFL:high.sps = PD:high.sps")

het2.sps <- felm(fu.sps.ind ~ SFL + PD + SFL:high.sps + 
                   PD:high.sps + bl.sps.ind + high.sps|0|0, data = full)
het2.sps2 <- felm(fu.sps.ind ~ SFL + PD + SFL:low.sps + 
                    PD:low.sps + bl.sps.ind + high.sps|0|0, data = full)
summary(het2.sps, robust = TRUE)
linearHypothesis(het2.sps, "SFL = PD")
linearHypothesis(het2.sps, "SFL:high.sps = PD:high.sps")

het2.e <- felm(fu.e.ind ~ SFL + PD + SFL:high.sps + 
                 PD:high.sps + bl.e.ind + high.sps|0|0, data = full)
het2.e2 <- felm(fu.e.ind ~ SFL + PD + SFL:low.sps + 
                  PD:low.sps + bl.e.ind + high.sps|0|0, data = full)
summary(het2.e, robust = TRUE)
linearHypothesis(het2.e, "SFL = PD")
linearHypothesis(het2.e, "SFL:high.sps = PD:high.sps")

coefs <- c(het2.pps$coefficients[2], het2.pps$coefficients[3],
           het2.pps2$coefficients[2], het2.pps2$coefficients[3],
           het2.sps$coefficients[2], het2.sps$coefficients[3],
           het2.sps2$coefficients[2], het2.sps2$coefficients[3],
           het2.e$coefficients[2], het2.e$coefficients[3],
           het2.e2$coefficients[2], het2.e2$coefficients[3])
SEs <- c(het2.pps$rse[2], het2.pps$rse[3],
         het2.pps2$rse[2], het2.pps2$rse[3],
         het2.sps$rse[2], het2.sps$rse[3],
         het2.sps2$rse[2], het2.sps2$rse[3],
         het2.e$rse[2], het2.e$rse[3],
         het2.e2$rse[2], het2.e2$rse[3])

het.sps.pdata <- as.data.frame(cbind(names, tmt, group, outcome, coefs, SEs))
het.sps.pdata$coefs <- as.numeric(het.sps.pdata$coefs)
het.sps.pdata$SEs <- as.numeric(het.sps.pdata$SEs)
het.sps.pdata <- mutate(het.sps.pdata, CI.High = coefs + 1.96*SEs,
                        CI.Low = coefs - 1.96*SEs)

#####Heterogeneity by above median econ index####
median.e <- median(full$bl.e.ind, na.rm=T)
full <- mutate(full, high.e = if_else(bl.e.ind>=median.e,1,0),
               low.e = if_else(bl.e.ind<median.e,1,0))

het3.pps <- felm(fu.pps.ind ~ SFL + PD + SFL:high.e + 
                   PD:high.e + bl.pps.ind + high.e|0|0, data = full)
het3.pps2 <- felm(fu.pps.ind ~ SFL + PD + SFL:low.e + 
                    PD:low.e + bl.pps.ind + high.e|0|0, data = full)
summary(het3.pps, robust = TRUE)
linearHypothesis(het3.pps, "SFL = PD")
linearHypothesis(het3.pps, "SFL:high.e = PD:high.e")

het3.sps <- felm(fu.sps.ind ~ SFL + PD + SFL:high.e + 
                   PD:high.e + bl.sps.ind + high.e|0|0, data = full)
het3.sps2 <- felm(fu.sps.ind ~ SFL + PD + SFL:low.e + 
                    PD:low.e + bl.sps.ind + high.e|0|0, data = full)
summary(het3.sps, robust = TRUE)
linearHypothesis(het3.sps, "SFL = PD")
linearHypothesis(het3.sps, "SFL:high.e = PD:high.e")

het3.e <- felm(fu.e.ind ~ SFL + PD + SFL:high.e + 
                 PD:high.e + bl.e.ind + high.e|0|0, data = full)
het3.e2 <- felm(fu.e.ind ~ SFL + PD + SFL:low.e + 
                  PD:low.e + bl.e.ind + high.e|0|0, data = full)
summary(het3.e, robust = TRUE)
linearHypothesis(het3.e, "SFL = PD")
linearHypothesis(het3.e, "SFL:high.e = PD:high.e")

coefs <- c(het3.pps$coefficients[2], het3.pps$coefficients[3],
           het3.pps2$coefficients[2], het3.pps2$coefficients[3],
           het3.sps$coefficients[2], het3.sps$coefficients[3],
           het3.sps2$coefficients[2], het3.sps2$coefficients[3],
           het3.e$coefficients[2], het3.e$coefficients[3],
           het3.e2$coefficients[2], het3.e2$coefficients[3])
SEs <- c(het3.pps$rse[2], het3.pps$rse[3],
         het3.pps2$rse[2], het3.pps2$rse[3],
         het3.sps$rse[2], het3.sps$rse[3],
         het3.sps2$rse[2], het3.sps2$rse[3],
         het3.e$rse[2], het3.e$rse[3],
         het3.e2$rse[2], het3.e2$rse[3])

het.e.pdata <- as.data.frame(cbind(names, tmt, group, outcome, coefs, SEs))
het.e.pdata$coefs <- as.numeric(het.e.pdata$coefs)
het.e.pdata$SEs <- as.numeric(het.e.pdata$SEs)
het.e.pdata <- mutate(het.e.pdata, CI.High = coefs + 1.96*SEs,
                      CI.Low = coefs - 1.96*SEs)

#####Indices of all pre-specified variables within each family####
#Primary psychosocial
pps.reg <- felm(fu.pps.ind ~ SFL + PD + bl.pps.ind |0|0, data = full)
summary(pps.reg, robust = TRUE)
linearHypothesis(pps.reg, "SFL = PD")

pps.reg2 <- felm(fu.pps.ind ~ SFL + PD |0|0, data = full)
summary(pps.reg2, robust = TRUE)
linearHypothesis(pps.reg2, "SFL = PD")

#Secondary psychosocial
sps.reg <- felm(fu.sps.ind ~ SFL + PD + bl.sps.ind|0|0, data = full)
summary(sps.reg, robust = TRUE)
linearHypothesis(sps.reg, "SFL = PD")

sps.reg2 <- felm(fu.sps.ind ~ SFL + PD |0|0, data = full)
summary(sps.reg2, robust = TRUE)
linearHypothesis(sps.reg2, "SFL = PD")

#Economic
e.reg <- felm(fu.e.ind ~ SFL + PD + bl.e.ind|0|0, data = full)
summary(e.reg, robust = TRUE)
linearHypothesis(e.reg, "SFL = PD")

e.reg2 <- felm(fu.e.ind ~ SFL + PD |0|0, data = full)
summary(e.reg2, robust = TRUE)
linearHypothesis(e.reg2, "SFL = PD")

#####Potential spillovers####
c.bytmt <- mutate(full, counter = 1) %>% group_by(Coop, workshop_name) %>% summarize(n = sum(counter))

c.sfl <- filter(c.bytmt, workshop_name=="Storytelling for Leadership" & !is.na(Coop)) %>% select(-workshop_name) %>%
  rename(c.SFL = n)
c.pd <- filter(c.bytmt, workshop_name=="Professional Development" & !is.na(Coop)) %>% select(-workshop_name) %>%
  rename(c.pd = n)

full <- left_join(full, c.sfl,by = "Coop")
full$c.SFL[which(is.na(full$c.SFL))] <- 0
full <- left_join(full, c.pd, by = "Coop")
full$c.pd[which(is.na(full$c.pd))] <- 0

median(full$c.SFL)
median(full$c.pd)

full <- mutate(full, c.sfl.high = if_else(c.SFL>3,1,0), c.pd.high = if_else(c.pd>3,1,0))

#Run everything only on the control group
c.pps.reg <- felm(fu.pps.ind ~ bl.pps.ind + c.SFL + c.pd|0|0, data = full[which(full$workshop_name=="Control"),])
summary(c.pps.reg, robust = TRUE)

c.sps.reg <- felm(fu.sps.ind ~ bl.sps.ind + c.SFL + c.pd|0|0, data = full[which(full$workshop_name=="Control"),])
summary(c.sps.reg, robust = TRUE)

c.e.reg <- felm(fu.e.ind ~ bl.e.ind + c.SFL + c.pd|0|0, data = full[which(full$workshop_name=="Control"),])
summary(c.e.reg, robust = TRUE)

#With above/below median dummy
c2.pps.reg <- felm(fu.pps.ind ~ bl.pps.ind + c.sfl.high + c.pd.high|0|0, data = full[which(full$workshop_name=="Control"),])
summary(c2.pps.reg, robust = TRUE)

c2.sps.reg <- felm(fu.sps.ind ~ bl.sps.ind + c.sfl.high + c.pd.high|0|0, data = full[which(full$workshop_name=="Control"),])
summary(c2.sps.reg, robust = TRUE)

c2.e.reg <- felm(fu.e.ind ~ bl.e.ind + c.sfl.high + c.pd.high|0|0, data = full[which(full$workshop_name=="Control"),])
summary(c2.e.reg, robust = TRUE)

#####Individual regressions for outcomes in the communication index####
# Create a separate dataset for the outcomes in this index
id.info <- select(full, resp_id, SFL, PD)
bl.com <- select(bl, resp_id, bl.suggestions = sei_6, bl.raise.ideas = sei_7a, 
                 bl.speak.up = sei_8a, bl.understand.suggestions = sei_16a, 
                 bl.understand.feedback = sei_16b, bl.understand.requests = sei_16c)
fu.com <- select(fu, resp_id, fu.suggestions = sei_6, fu.raise.ideas = sei_7a, 
                 fu.speak.up = sei_8a, fu.understand.suggestions = sei_16a, 
                 fu.understand.feedback = sei_16b, fu.understand.requests = sei_16c)
com <- left_join(id.info, bl.com)

com <- left_join(com, fu.com)

creg1 <- felm(fu.suggestions ~ SFL + PD + bl.suggestions|0|0, data = com)
summary(creg1, robust = TRUE)
linearHypothesis(creg1, "SFL = PD")

creg2 <- felm(fu.raise.ideas ~ SFL + PD + bl.raise.ideas|0|0, data = com)
summary(creg2, robust = TRUE)
linearHypothesis(creg2, "SFL = PD")

creg3 <- felm(fu.speak.up ~ SFL + PD + bl.speak.up|0|0, data = com)
summary(creg3, robust = TRUE)
linearHypothesis(creg3, "SFL = PD")

creg4 <- felm(fu.understand.suggestions ~ SFL + PD + bl.understand.suggestions|0|0, data = com)
summary(creg4, robust = TRUE)
linearHypothesis(creg4, "SFL = PD")

creg5 <- felm(fu.understand.feedback ~ SFL + PD + bl.understand.feedback|0|0, data = com)
summary(creg5, robust = TRUE)
linearHypothesis(creg5, "SFL = PD")

creg6 <- felm(fu.understand.requests ~ SFL + PD + bl.understand.requests|0|0, data = com)
summary(creg6, robust = TRUE)
linearHypothesis(creg6, "SFL = PD")

##########################################################################################
#Control the false discovery rate within families of outcomes
##########################################################################################
##### Primary psychosocial ####
SFL.pps <- as.data.frame(cbind(c("Self-Advocacy Score","Problem Solving Score","No. Leadership Positions", "Connectedness to Home/Village", 
                                 "Connectedness to Co-op", "Communication Skills Score", "Any Leadership Position"), 
                               c(sa.reg.ew$coefficients[2], ps.reg.ew$coefficients[2], nlead.reg$coefficients[2], conhv.reg.ew$coefficients[2], 
                                 con.reg.ew$coefficients[2],com.reg.ew$coefficients[2], lead.ind.reg$coefficients[2]),
                               c(sa.reg.ew$rpval[2], ps.reg.ew$rpval[2], nlead.reg$rpval[2], conhv.reg.ew$rpval[2], con.reg.ew$rpval[2],
                                 com.reg.ew$rpval[2], lead.ind.reg$rpval[2]))) 
SFL.pps$V2 <- as.numeric(SFL.pps$V2)
SFL.pps$V3 <- as.numeric(SFL.pps$V3)
colnames(SFL.pps) <- c("Outcome","coef","pval")
SFL.pps <- gen.qvals(SFL.pps)

PD.pps <- as.data.frame(cbind(c("Self-Advocacy Score","Problem Solving Score","No. Leadership Positions", "Connectedness to Home/Village", 
                                "Connectedness to Co-op", "Communication Skills Score", "Any Leadership Position"), 
                              c(sa.reg.ew$coefficients[3], ps.reg.ew$coefficients[3], nlead.reg$coefficients[3], conhv.reg.ew$coefficients[3], 
                                con.reg.ew$coefficients[3],com.reg.ew$coefficients[3], lead.ind.reg$coefficients[3]),
                              c(sa.reg.ew$rpval[3], ps.reg.ew$rpval[3], nlead.reg$rpval[3], conhv.reg.ew$rpval[3], con.reg.ew$rpval[3],
                                com.reg.ew$rpval[3], lead.ind.reg$rpval[3]))) 
PD.pps$V2 <- as.numeric(PD.pps$V2)
PD.pps$V3 <- as.numeric(PD.pps$V3)
colnames(PD.pps) <- c("Outcome","coef","pval")

PD.pps <- gen.qvals(PD.pps)

##### Secondary psychosocial ####
SFL.sps <- as.data.frame(cbind(c("SSS - Current","SSS - Past","SSS - Future","Connectedness to Workshop Group", "Self-Value Score",
                                 "Resilience Score","ID as Leader Score","Has Role Model", "Aspires to be Role Model",
                                 "Time to be Role Model","Has a Goal","Achieved Last Goal"),
                               c(ss.reg$coefficients[2], ssp.reg$coefficients[2], ssf.reg$coefficients[2], tg.reg$coefficients[2],
                                 val.reg.ew$coefficients[2], res.reg.ew$coefficients[2], lead.reg$coefficients[2], rm.reg$coefficients[2],
                                 rmasp.reg$coefficients[2], rmtime.reg$coefficients[2], goal.reg$coefficients[2], com.goal.reg$coefficients[2]),
                               c(ss.reg$rpval[2], ssp.reg$rpval[2], ssf.reg$rpval[2], tg.reg$rpval[2], val.reg.ew$rpval[2], res.reg.ew$rpval[2], 
                                 lead.reg$rpval[2], rm.reg$rpval[2], rmasp.reg$rpval[2], rmtime.reg$rpval[2], goal.reg$rpval[2], 
                                 com.goal.reg$rpval[2])))

SFL.sps$V2 <- as.numeric(SFL.sps$V2)
SFL.sps$V3 <- as.numeric(SFL.sps$V3)
colnames(SFL.sps) <- c("Outcome","coef","pval")
SFL.sps <- gen.qvals(SFL.sps)

PD.sps <- as.data.frame(cbind(c("SSS - Current","SSS - Past","SSS - Future","Connectedness to Workshop Group", "Self-Value Score",
                                "Resilience Score","ID as Leader Score","Has Role Model", "Aspires to be Role Model",
                                "Time to be Role Model","Has a Goal","Achieved Last Goal"),
                              c(ss.reg$coefficients[3], ssp.reg$coefficients[3], ssf.reg$coefficients[3], tg.reg$coefficients[3],
                                val.reg.ew$coefficients[3], res.reg.ew$coefficients[3], lead.reg$coefficients[3], rm.reg$coefficients[3],
                                rmasp.reg$coefficients[3], rmtime.reg$coefficients[3], goal.reg$coefficients[3], com.goal.reg$coefficients[3]),
                              c(ss.reg$rpval[3], ssp.reg$rpval[3], ssf.reg$rpval[3], tg.reg$rpval[3], val.reg.ew$rpval[3], res.reg.ew$rpval[3], 
                                lead.reg$rpval[3], rm.reg$rpval[3], rmasp.reg$rpval[3], rmtime.reg$rpval[3], goal.reg$rpval[3], 
                                com.goal.reg$rpval[3])))

PD.sps$V2 <- as.numeric(PD.sps$V2)
PD.sps$V3 <- as.numeric(PD.sps$V3)
colnames(PD.sps) <- c("Outcome","coef","pval")
PD.sps <- gen.qvals(PD.sps)

##### Economic ####
SFL.e <- as.data.frame(cbind(c("Marginal Utility of Expenditure","Log(Consumption Expenditures per Capita)","Earn Any Income","Monthly Cash Earnings (RWF)"),
                             c(mue.reg$coefficients[2], exp.reg$coefficients[2], work.reg$coefficients[2], inc$coefficients[2]),
                             c(mue.reg$rpval[2], exp.reg$rpval[2], work.reg$rpval[2], inc$rpval[2])))

SFL.e$V2 <- as.numeric(SFL.e$V2)
SFL.e$V3 <- as.numeric(SFL.e$V3)
colnames(SFL.e) <- c("Outcome","coef","pval")
SFL.e <- gen.qvals(SFL.e)

PD.e <- as.data.frame(cbind(c("Marginal Utility of Expenditure","Log(Consumption Expenditures per Capita)","Earn Any Income","Monthly Cash Earnings (RWF)"),
                            c(mue.reg$coefficients[3], exp.reg$coefficients[3], work.reg$coefficients[3], inc$coefficients[3]),
                            c(mue.reg$rpval[3], exp.reg$rpval[3], work.reg$rpval[3], inc$rpval[3])))

PD.e$V2 <- as.numeric(PD.e$V2)
PD.e$V3 <- as.numeric(PD.e$V3)
colnames(PD.e) <- c("Outcome","coef","pval")
PD.e <- gen.qvals(PD.e)

##########################################################################################
#Figures 
##########################################################################################
#####Figure 1 - Treatment effect plots for primary psycho-social outcomes#####
sa.sd <- sd(full$fu.saew[which(full$workshop_name=="Control")], na.rm=T)
sa.mean <- mean(full$fu.saew[which(full$workshop_name=="Control")], na.rm=T)
sa.pd <- as.data.frame(cbind(sa.reg.ew$coefficients, sa.reg.ew$rse))
sa.pd <- cbind(sa.pd, rownames(sa.pd))
colnames(sa.pd) <- c("ATE","RSE","Tmt")
sa.pd <- mutate(sa.pd, dep.var = "Self-Advocacy Score")

l.sd <- sd(full$fu.lead[which(full$workshop_name=="Control")], na.rm=T)
l.mean <- mean(full$fu.lead[which(full$workshop_name=="Control")], na.rm=T)
l.pd <- as.data.frame(cbind(nlead.reg$coefficients, nlead.reg$rse))
l.pd <- cbind(l.pd, rownames(l.pd))
colnames(l.pd) <- c("ATE","RSE","Tmt")
l.pd <- mutate(l.pd, dep.var = "No. Leadership Positions")

lind.sd <- sqrt((mean(full$fu.lead.ind[which(full$workshop_name=="Control")], na.rm=T))*(1-mean(full$fu.lead.ind[which(full$workshop_name=="Control")], na.rm=T)))
lind.mean <- mean(full$fu.lead.ind[which(full$workshop_name=="Control")], na.rm=T)
lind.pd <- as.data.frame(cbind(lead.ind.reg$coefficients, lead.ind.reg$rse))
lind.pd <- cbind(lind.pd, rownames(lind.pd))
colnames(lind.pd) <- c("ATE","RSE","Tmt")
lind.pd <- mutate(lind.pd, dep.var = "Any Leadership Position")

cc.sd <- sd(full$fu.conew[which(full$workshop_name=="Control")], na.rm=T)
cc.mean <- mean(full$fu.conew[which(full$workshop_name=="Control")], na.rm=T)
cc.pd <- as.data.frame(cbind(con.reg.ew$coefficients, con.reg.ew$rse))
cc.pd <- cbind(cc.pd, rownames(cc.pd))
colnames(cc.pd) <- c("ATE","RSE","Tmt")
cc.pd <- mutate(cc.pd, dep.var = "Connectedness to Co-op")

chv.sd <- sd(full$fu.conhvew[which(full$workshop_name=="Control")], na.rm=T)
chv.mean <- mean(full$fu.conhvew[which(full$workshop_name=="Control")], na.rm=T)
chv.pd <- as.data.frame(cbind(conhv.reg.ew$coefficients, conhv.reg.ew$rse))
chv.pd <- cbind(chv.pd, rownames(chv.pd))
colnames(chv.pd) <- c("ATE","RSE","Tmt")
chv.pd <- mutate(chv.pd, dep.var = "Connectedness to Home/Village")

ps.sd <- sd(full$fu.psew[which(full$workshop_name=="Control")], na.rm=T)
ps.mean <- mean(full$fu.psew[which(full$workshop_name=="Control")], na.rm=T)
ps.pd <- as.data.frame(cbind(ps.reg.ew$coefficients, ps.reg.ew$rse))
ps.pd <- cbind(ps.pd, rownames(ps.pd))
colnames(ps.pd) <- c("ATE","RSE","Tmt")
ps.pd <- mutate(ps.pd, dep.var = "Problem Solving Score")

com.sd <- sd(full$fu.comew[which(full$workshop_name=="Control")], na.rm=T)
com.mean <- mean(full$fu.comew[which(full$workshop_name=="Control")], na.rm=T)
com.pd <- as.data.frame(cbind(com.reg.ew$coefficients, com.reg.ew$rse))
com.pd <- cbind(com.pd, rownames(com.pd))
colnames(com.pd) <- c("ATE","RSE","Tmt")
com.pd <- mutate(com.pd, dep.var = "Communication Skills Score")

pps.sd <- sd(full$fu.pps.ind[which(full$workshop_name=="Control")], na.rm=T)
pps.mean <- mean(full$fu.pps.ind[which(full$workshop_name=="Control")], na.rm=T)
pps.pd <- as.data.frame(cbind(pps.reg$coefficients, pps.reg$rse))
pps.pd <- cbind(pps.pd, rownames(pps.pd))
colnames(pps.pd) <- c("ATE","RSE","Tmt")
pps.pd <- mutate(pps.pd, dep.var = "PPS Index")

pps.plot.data <- rbind(sa.pd, l.pd, lind.pd, cc.pd, chv.pd, ps.pd, com.pd, pps.pd) 
pps.plot.data <- filter(pps.plot.data, Tmt=="SFL"|Tmt=="PD") %>%
  mutate(SD = if_else(dep.var=="Self-Advocacy Score", sa.sd, 
                      if_else(dep.var=="Connectedness to Home/Village", chv.sd,
                              if_else(dep.var=="Problem Solving Score", ps.sd, 
                                      if_else(dep.var=="Communication Skills Score", com.sd,
                                              if_else(dep.var=="Connectedness to Co-op",cc.sd,
                                                      if_else(dep.var=="No. Leadership Positions",l.sd,
                                                              if_else(dep.var=="Any Leadership Position", lind.sd, pps.sd))))))))

#Merge in multiple inference-adjusted q-values and SEs
SFL.pps <- mutate(SFL.pps, Tmt = "SFL") %>% select("Outcome","qval","adj.SE","Tmt")
PD.pps <- mutate(PD.pps, Tmt = "PD") %>% select("Outcome","qval","adj.SE","Tmt")
pps.adj <- rbind(SFL.pps, PD.pps)
pps.plot.data <- left_join(pps.plot.data, pps.adj, by = c("dep.var" = "Outcome", "Tmt"))
pps.plot.data <- mutate(pps.plot.data, ATE.std = ATE/SD, RSE.std = RSE/SD, CI.high = ATE.std + 1.965*RSE.std, CI.low = ATE.std - 1.965*RSE.std,
                        adj.RSE.std = adj.SE/SD, adj.CI.high = ATE.std + 1.965*adj.RSE.std, adj.CI.low = ATE.std - 1.965*adj.RSE.std)   
pps.plot.data$dep.var[which(pps.plot.data$dep.var=="PPS Index")] <- " PPS Index"
pps.plot.data <- mutate(pps.plot.data, Treatment = if_else(Tmt=="PD","Soft Skills","Storytelling"))

pps.plot <- ggplot(data = pps.plot.data, aes(y = dep.var, x = ATE.std, group = Treatment,
                                             fill = Treatment, color = Treatment)) + 
  geom_point(position = position_dodge(0.5)) + 
  geom_errorbar(aes(xmin = CI.low, xmax = CI.high), width = 0.4, position = position_dodge(0.5)) + 
  geom_errorbar(aes(xmin = adj.CI.low, xmax = adj.CI.high), width = 0.2, position = position_dodge(0.5)) + 
  geom_vline(xintercept = 0, color = "black") + 
  labs(x = "Treatment Effect (Standard Deviations)", 
       y = "Outcome") + 
  scale_fill_grey(start = 0.5, end = 0.05) + scale_color_grey(start = 0.5, end = 0.05) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) 

ggsave(filename="../Figures/PPS.eps", plot=pps.plot, device="eps", height=2.12, width=5, units="in", dpi = 400)

#####Figure 2 - Treatment effect plots for secondary psycho-social outcomes#####
pss.sd <- sd(full$fup.ss[which(full$workshop_name=="Control")], na.rm=T)
pss.mean <- mean(full$fup.ss[which(full$workshop_name=="Control")], na.rm=T)
pss.pd <- as.data.frame(cbind(ssp.reg$coefficients, ssp.reg$rse))
pss.pd <- cbind(pss.pd, rownames(pss.pd))
colnames(pss.pd) <- c("ATE","RSE","Tmt")
pss.pd <- mutate(pss.pd, dep.var = "SSS - Past")

ss.sd <- sd(full$fu.ss[which(full$workshop_name=="Control")], na.rm=T)
ss.mean <- mean(full$fu.ss[which(full$workshop_name=="Control")], na.rm=T)
ss.pd <- as.data.frame(cbind(ss.reg$coefficients, ss.reg$rse))
ss.pd <- cbind(ss.pd, rownames(ss.pd))
colnames(ss.pd) <- c("ATE","RSE","Tmt")
ss.pd <- mutate(ss.pd, dep.var = "SSS - Current")

fss.sd <- sd(full$fuf.ss[which(full$workshop_name=="Control")], na.rm=T)
fss.mean <- mean(full$fuf.ss[which(full$workshop_name=="Control")], na.rm=T)
fss.pd <- as.data.frame(cbind(ssf.reg$coefficients, ssf.reg$rse))
fss.pd <- cbind(fss.pd, rownames(fss.pd))
colnames(fss.pd) <- c("ATE","RSE","Tmt")
fss.pd <- mutate(fss.pd, dep.var = "SSS - Future")

conw.sd <- sd(full$fu.tg[which(full$workshop_name=="Control")],na.rm=T)
conw.mean <- mean(full$fu.tg[which(full$workshop_name=="Control")],na.rm=T)
conw.pd <- as.data.frame(cbind(tg.reg$coefficients, tg.reg$rse))
conw.pd <- cbind(conw.pd, rownames(conw.pd))
colnames(conw.pd) <- c("ATE","RSE","Tmt")
conw.pd <- mutate(conw.pd, dep.var = "Connectedness to Workshop Group")

v.sd <- sd(full$fu.vew[which(full$workshop_name=="Control")],na.rm=T)
v.mean <- mean(full$fu.vew[which(full$workshop_name=="Control")],na.rm=T)
v.pd <- as.data.frame(cbind(val.reg.ew$coefficients, val.reg.ew$rse))
v.pd <- cbind(v.pd, rownames(v.pd))
colnames(v.pd) <- c("ATE","RSE","Tmt")
v.pd <- mutate(v.pd, dep.var = "Self-Value Score")

r.sd <- sd(full$fu.rew[which(full$workshop_name=="Control")], na.rm=T)
r.mean <- mean(full$fu.rew[which(full$workshop_name=="Control")], na.rm=T)
r.pd <- as.data.frame(cbind(res.reg.ew$coefficients, res.reg.ew$rse))
r.pd <- cbind(r.pd, rownames(r.pd))
colnames(r.pd) <- c("ATE","RSE","Tmt")
r.pd <- mutate(r.pd, dep.var = "Resilience Score")

idl.sd <- sd(full$fu.lead[which(full$workshop_name=="Control")], na.rm=T)
idl.mean <- mean(full$fu.lead[which(full$workshop_name=="Control")], na.rm=T)
idl.pd <- as.data.frame(cbind(lead.reg$coefficients, lead.reg$rse))
idl.pd <- cbind(idl.pd, rownames(idl.pd))
colnames(idl.pd) <- c("ATE","RSE","Tmt")
idl.pd <- mutate(idl.pd, dep.var = "ID as Leader Score")

rm.sd <- sqrt(mean(full$fu.rm[which(full$workshop_name=="Control")], na.rm=T)*(1-mean(full$fu.rm[which(full$workshop_name=="Control")], na.rm=T)))
rm.mean <- mean(full$fu.rm[which(full$workshop_name=="Control")], na.rm=T)
rm.pd <- as.data.frame(cbind(rm.reg$coefficients, rm.reg$rse))
rm.pd <- cbind(rm.pd, rownames(rm.pd))
colnames(rm.pd) <- c("ATE","RSE","Tmt")
rm.pd <- mutate(rm.pd, dep.var = "Has Role Model")

rma.sd <- sqrt(mean(full$fu.rmasp[which(full$workshop_name=="Control")], na.rm=T)*(1-mean(full$fu.rmasp[which(full$workshop_name=="Control")], na.rm=T)))
rma.mean <- mean(full$fu.rmasp[which(full$workshop_name=="Control")], na.rm=T)
rma.pd <- as.data.frame(cbind(rmasp.reg$coefficients, rmasp.reg$rse))
rma.pd <- cbind(rma.pd, rownames(rma.pd))
colnames(rma.pd) <- c("ATE","RSE","Tmt")
rma.pd <- mutate(rma.pd, dep.var = "Aspires to be Role Model")

rmt.sd <- sd(full$fu.rmtime[which(full$workshop_name=="Control")], na.rm=T)
rmt.mean <- mean(full$fu.rmtime[which(full$workshop_name=="Control")], na.rm=T)
rmt.pd <- as.data.frame(cbind(rmtime.reg$coefficients, rmtime.reg$rse))
rmt.pd <- cbind(rmt.pd, rownames(rmt.pd))
colnames(rmt.pd) <- c("ATE","RSE","Tmt")
rmt.pd <- mutate(rmt.pd, dep.var = "Time to be Role Model")

g.sd <- sqrt(mean(full$fu.goal[which(full$workshop_name=="Control")], na.rm=T)*(1-mean(full$fu.goal[which(full$workshop_name=="Control")], na.rm=T)))
g.mean <- mean(full$fu.goal[which(full$workshop_name=="Control")], na.rm=T)
g.pd <- as.data.frame(cbind(goal.reg$coefficients, goal.reg$rse))
g.pd <- cbind(g.pd, rownames(g.pd))
colnames(g.pd) <- c("ATE","RSE","Tmt")
g.pd <- mutate(g.pd, dep.var = "Has a Goal")

gc.sd <- sqrt(mean(full$com.goal[which(full$workshop_name=="Control")], na.rm=T)*(1-mean(full$com.goal[which(full$workshop_name=="Control")], na.rm=T)))
gc.mean <- mean(full$com.goal[which(full$workshop_name=="Control")], na.rm=T)
gc.pd <- as.data.frame(cbind(com.goal.reg$coefficients, com.goal.reg$rse))
gc.pd <- cbind(gc.pd, rownames(gc.pd))
colnames(gc.pd) <- c("ATE","RSE","Tmt")
gc.pd <- mutate(gc.pd, dep.var = "Achieved Last Goal")

sps.sd <- sd(full$fu.sps.ind[which(full$workshop_name=="Control")], na.rm=T)
sps.mean <- mean(full$fu.sps.ind[which(full$workshop_name=="Control")], na.rm=T)
sps.pd <- as.data.frame(cbind(sps.reg$coefficients, sps.reg$rse))
sps.pd <- cbind(sps.pd, rownames(sps.pd))
colnames(sps.pd) <- c("ATE","RSE","Tmt")
sps.pd <- mutate(sps.pd, dep.var = "SPS Index")

sps.plot.data <- rbind(pss.pd, ss.pd, fss.pd, conw.pd, v.pd, r.pd, idl.pd, rm.pd, rma.pd, rmt.pd, 
                       g.pd, gc.pd, sps.pd) 
sps.plot.data <- filter(sps.plot.data, Tmt=="SFL"|Tmt=="PD") %>%
  mutate(SD = if_else(dep.var=="SSS - Past", pss.sd, 
                      if_else(dep.var=="SSS - Current", ss.sd,
                              if_else(dep.var=="SSS - Future", fss.sd, 
                                      if_else(dep.var=="Connectedness to Workshop Group", conw.sd,
                                              if_else(dep.var=="Self-Value Score",v.sd,
                                                      if_else(dep.var=="Resilience Score",r.sd,
                                                              if_else(dep.var=="ID as Leader Score", idl.sd,
                                                                      if_else(dep.var=="Has Role Model", rm.sd,
                                                                              if_else(dep.var=="Aspires to be Role Model", rma.sd,
                                                                                      if_else(dep.var=="Time to be Role Model", rmt.sd,
                                                                                              if_else(dep.var=="Has a Goal", g.sd,
                                                                                                      if_else(dep.var=="Achieved Last Goal", gc.sd,
                                                                                                              sps.sd)))))))))))))
#Merge in multiple inference-adjusted q-values and SEs
SFL.sps <- mutate(SFL.sps, Tmt = "SFL") %>% select("Outcome","qval","adj.SE","Tmt")
PD.sps <- mutate(PD.sps, Tmt = "PD") %>% select("Outcome","qval","adj.SE","Tmt")
sps.adj <- rbind(SFL.sps, PD.sps)
sps.plot.data <- left_join(sps.plot.data, sps.adj, by = c("dep.var" = "Outcome", "Tmt"))
sps.plot.data <- mutate(sps.plot.data, ATE.std = ATE/SD, RSE.std = RSE/SD, CI.high = ATE.std + 1.965*RSE.std, CI.low = ATE.std - 1.965*RSE.std,
                        adj.RSE.std = adj.SE/SD, adj.CI.high = ATE.std + 1.965*adj.RSE.std, adj.CI.low = ATE.std - 1.965*adj.RSE.std)   

sps.plot.data$dep.var <- factor(sps.plot.data$dep.var, 
                                levels = c("SPS Index","Achieved Last Goal","Has a Goal",
                                           "Time to be Role Model","Aspires to be Role Model",
                                           "Has Role Model","ID as Leader Score","Resilience Score", 
                                           "Self-Value Score", "Connectedness to Workshop Group", 
                                           "SSS - Future","SSS - Past","SSS - Current"))
sps.plot.data <- mutate(sps.plot.data, Treatment = if_else(Tmt=="PD","Soft Skills","Storytelling"))

sps.plot.1 <- ggplot(data = sps.plot.data, 
                     aes(y = dep.var, x = ATE.std, group = Treatment,
                         fill = Treatment, color = Treatment)) + 
  geom_point(position = position_dodge(0.5)) + coord_cartesian(xlim=c(-0.5,1.75)) + 
  geom_errorbar(aes(xmin = CI.low, xmax = CI.high), width = 0.3, position = position_dodge(0.5)) + 
  geom_errorbar(aes(xmin = adj.CI.low, xmax = adj.CI.high), width = 0.15, position = position_dodge(0.5)) + 
  geom_vline(xintercept = 0, color = "black") + 
  labs(x = "Treatment Effect (Standard Deviations)", 
       y = "Outcome") + 
  scale_fill_grey(start = 0.5, end = 0.05) + scale_color_grey(start = 0.5, end = 0.05) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5))

ggsave(filename="../Figures/SPS.eps", plot=sps.plot.1, device="eps", height=3, width=5, units="in", dpi=400)

#####Figure 3 - Treatment effect plots for economic outcomes#####
inc.sd <- sd(full$fu.income[which(full$workshop_name=="Control")], na.rm=T)
inc.mean <- mean(full$fu.income[which(full$workshop_name=="Control")], na.rm=T)
inc.pd <- as.data.frame(cbind(inc$coefficients, inc$rse))
inc.pd <- cbind(inc.pd, rownames(inc.pd))
colnames(inc.pd) <- c("ATE","RSE","Tmt")
inc.pd <- mutate(inc.pd, dep.var = "Monthly Cash Earnings (RWF)")

inci.sd <- sqrt(mean(full$fu.work_lastmonth[which(full$workshop_name=="Control")], na.rm=T)*(1-mean(full$fu.work_lastmonth[which(full$workshop_name=="Control")], na.rm=T)))
inci.mean <- mean(full$fu.work_lastmonth[which(full$workshop_name=="Control")], na.rm=T)
inci.pd <- as.data.frame(cbind(work.reg$coefficients, work.reg$rse))
inci.pd <- cbind(inci.pd, rownames(inci.pd))
colnames(inci.pd) <- c("ATE","RSE","Tmt")
inci.pd <- mutate(inci.pd, dep.var = "Earn Any Income")

mue.sd <- sd(full$fu.mue[which(full$workshop_name=="Control")], na.rm=T)
mue.mean <- mean(full$fu.mue[which(full$workshop_name=="Control")], na.rm=T)
mue.pd <- as.data.frame(cbind(mue.reg$coefficients, mue.reg$rse))
mue.pd <- cbind(mue.pd, rownames(mue.pd))
colnames(mue.pd) <- c("ATE","RSE","Tmt")
mue.pd <- mutate(mue.pd, dep.var = "Marginal Utility of Expenditure")

exp.sd <- sd(full$fu.log.exp.pc[which(full$workshop_name=="Control")], na.rm=T)
exp.mean <- mean(full$fu.log.exp.pc[which(full$workshop_name=="Control")], na.rm=T)
exp.pd <- as.data.frame(cbind(exp.reg$coefficients, exp.reg$rse))
exp.pd <- cbind(exp.pd, rownames(exp.pd))
colnames(exp.pd) <- c("ATE","RSE","Tmt")
exp.pd <- mutate(exp.pd, dep.var = "Log(Consumption Expenditures per Capita)")

e.sd <- sd(full$fu.e.ind[which(full$workshop_name=="Control")], na.rm=T)
e.mean <- mean(full$fu.e.ind[which(full$workshop_name=="Control")], na.rm=T)
e.pd <- as.data.frame(cbind(e.reg$coefficients, e.reg$rse))
e.pd <- cbind(e.pd, rownames(e.pd))
colnames(e.pd) <- c("ATE","RSE","Tmt")
e.pd <- mutate(e.pd, dep.var = "Econ Index")

e.plot.data <- rbind(inc.pd, inci.pd, mue.pd, exp.pd, e.pd) 
e.plot.data <- filter(e.plot.data, Tmt=="SFL"|Tmt=="PD") %>%
  mutate(SD = if_else(dep.var=="Monthly Cash Earnings (RWF)", inc.sd, 
                      if_else(dep.var=="Earn Any Income", inci.sd,
                              if_else(dep.var=="Marginal Utility of Expenditure", mue.sd,
                                      if_else(dep.var=="Log(Consumption Expenditures per Capita)",exp.sd, e.sd)))))

#Merge in multiple inference-adjusted q-values and SEs
SFL.e <- mutate(SFL.e, Tmt = "SFL") %>% select("Outcome","qval","adj.SE","Tmt")
PD.e <- mutate(PD.e, Tmt = "PD") %>% select("Outcome","qval","adj.SE","Tmt")
e.adj <- rbind(SFL.e, PD.e)
e.plot.data <- left_join(e.plot.data, e.adj, by = c("dep.var" = "Outcome", "Tmt"))
e.plot.data <- mutate(e.plot.data, ATE.std = ATE/SD, RSE.std = RSE/SD, CI.high = ATE.std + 1.965*RSE.std, CI.low = ATE.std - 1.965*RSE.std,
                      adj.RSE.std = adj.SE/SD, adj.CI.high = ATE.std + 1.965*adj.RSE.std, adj.CI.low = ATE.std - 1.965*adj.RSE.std)  
e.plot.data$dep.var[which(e.plot.data$dep.var=="Econ Index")] <- " Econ Index"
e.plot.data$dep.var[which(e.plot.data$dep.var=="Earn Any Income")] <- "Earn Any Cash"
e.plot.data <- mutate(e.plot.data, Treatment = if_else(Tmt=="PD","Soft Skills","Storytelling"))

e.plot <- ggplot(data = e.plot.data, aes(y = dep.var, x = ATE.std, group = Treatment,
                                         fill = Treatment, color = Treatment)) + 
  geom_point(position = position_dodge(0.5)) + 
  geom_errorbar(aes(xmin = CI.low, xmax = CI.high), width = 0.3, position = position_dodge(0.5)) + 
  geom_errorbar(aes(xmin = adj.CI.low, xmax = adj.CI.high), width = 0.15, position = position_dodge(0.5)) + 
  geom_vline(xintercept = 0, color = "black") + 
  labs(x = "Treatment Effect (Standard Deviations)", 
       y = "Outcome") + 
  scale_fill_grey(start = 0.5, end = 0.05) + scale_color_grey(start = 0.5, end = 0.05) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) 

ggsave(filename="../Figures/Econ.eps", plot=e.plot, device="eps", height=2.12, width=5, units="in", dpi=400)

##########################################################################################
#Tables 
##########################################################################################
#####Table 1 - Balance####
covs <- select(full, workshop_name, age, hhsize, bl.saew, bl.nlead, bl.lead.ind, bl.conew, bl.conhvew, 
               bl.psew, bl.comew, bl.income, bl.work_lastmonth, bl.mue, bl.rew, bl.vew,
               bl.rm, bl.rmasp, bl.rmtime, bl.ss, blp.ss, blf.ss, bl.lead, bl.goal)

bal.sa <- felm(bl.saew ~ SFL + PD|0|0, data = full)
summary(bal.sa, robust = TRUE)

bal.ps <- felm(bl.psew ~ SFL + PD|0|0, data = full)
summary(bal.ps, robust = TRUE)

bal.nlead <- felm(bl.nlead ~ SFL + PD|0|0, data = full)
summary(bal.nlead, robust = TRUE)

bal.conhv <- felm(bl.conhvew ~ SFL + PD|0|0, data = full)
summary(bal.conhv, robust = TRUE)

bal.con <- felm(bl.conew ~ SFL + PD|0|0, data = full)
summary(bal.con, robust = TRUE)

bal.com <- felm(bl.comew ~ SFL + PD|0|0, data = full)
summary(bal.com, robust = TRUE)

bal.lead.ind <- felm(bl.lead.ind ~ SFL + PD|0|0, data = full)
summary(bal.lead.ind, robust = TRUE)

bal.ss <- felm(bl.ss ~ SFL + PD|0|0, data = full)
summary(bal.ss, robust = TRUE)

bal.pss <- felm(blp.ss ~ SFL + PD|0|0, data = full)
summary(bal.pss, robust = TRUE)

bal.fss <- felm(blf.ss ~ SFL + PD|0|0, data = full)
summary(bal.fss, robust = TRUE)

bal.v <- felm(bl.vew ~ SFL + PD|0|0, data = full)
summary(bal.v, robust = TRUE)

bal.r <- felm(bl.rew ~ SFL + PD|0|0, data = full)
summary(bal.r, robust = TRUE)

bal.lead <- felm(bl.lead ~ SFL + PD|0|0, data = full)
summary(bal.lead, robust = TRUE)

bal.tg <- felm(bl.tg ~ SFL + PD|0|0, data = full)
summary(bal.tg, robust = TRUE)

bal.rm <- felm(bl.rm ~ SFL + PD|0|0, data = full)
summary(bal.rm, robust = TRUE)

bal.rmasp <- felm(bl.rmasp ~ SFL + PD|0|0, data = full)
summary(bal.rmasp, robust = TRUE)

bal.rmtime <- felm(bl.rmtime ~ SFL + PD|0|0, data = full)
summary(bal.rmtime, robust = TRUE)

bal.goal <- felm(bl.goal ~ SFL + PD|0|0, data = full)
summary(bal.goal, robust = TRUE)

bal.age <- felm(age ~ SFL + PD|0|0, data = full)
summary(bal.age, robust = TRUE)

bal.hhsize <- felm(hhsize ~ SFL + PD|0|0, data = full)
summary(bal.hhsize, robust = TRUE)

bal.inc <- felm(bl.income ~ SFL + PD|0|0, data = full)
summary(bal.inc, robust = TRUE)

bal.mue <- felm(bl.mue ~ SFL + PD|0|0, data = full)
summary(bal.mue, robust = TRUE)

bal.work <- felm(bl.work_lastmonth ~ SFL + PD|0|0, data = full)
summary(bal.work, robust = TRUE)

bal.cons <- felm(bl.log.exp.pc ~ SFL + PD|0|0, data = full)
summary(bal.cons, robust = TRUE)

#Make a table with means and standard errors
balance.means <- select(full, workshop_name, age, hhsize, bl.saew, bl.nlead, bl.lead.ind, bl.conew, 
                        bl.conhvew, bl.psew, bl.comew, bl.tg, bl.income, bl.work_lastmonth, bl.mue, bl.rew, bl.vew,
                        bl.rm, bl.rmasp, bl.rmtime, bl.ss, blp.ss, blf.ss, bl.lead, bl.goal, bl.log.exp.pc) %>%
  group_by(workshop_name) %>% summarize(mean.age = mean(age, na.rm=T), 
                                        mean.hhsize = mean(hhsize, na.rm=T), 
                                        mean.sa = mean(bl.saew, na.rm=T),
                                        mean.ps = mean(bl.psew, na.rm=T),
                                        mean.nlead = mean(bl.nlead, na.rm=T),
                                        mean.conhv = mean(bl.conhvew, na.rm=T),
                                        mean.con = mean(bl.conew, na.rm=T),
                                        mean.com = mean(bl.comew, na.rm=T),
                                        mean.lead.ind = mean(bl.lead.ind, na.rm=T),
                                        mean.ss = mean(bl.ss, na.rm=T),
                                        mean.pss = mean(blp.ss, na.rm=T),
                                        mean.fss = mean(blf.ss, na.rm=T),
                                        mean.v = mean(bl.vew, na.rm=T),
                                        mean.r = mean(bl.rew, na.rm=T),
                                        mean.lead = mean(bl.lead, na.rm=T),
                                        mean.tg = mean(bl.tg, na.rm=T),
                                        mean.rm = mean(bl.rm, na.rm=T),
                                        mean.rmasp = mean(bl.rmasp, na.rm=T),
                                        mean.rmtime = mean(bl.rmtime, na.rm=T),
                                        mean.goal = mean(bl.goal, na.rm=T),
                                        mean.income = mean(bl.income, na.rm=T), 
                                        mean.work = mean(bl.work_lastmonth, na.rm=T),
                                        mean.mue = mean(bl.mue, na.rm=T),
                                        mean.cons = mean(bl.log.exp.pc, na.rm=T))

balance.sds <- select(full, workshop_name, age, hhsize, bl.saew, bl.nlead, bl.lead.ind, bl.conew, 
                      bl.conhvew, bl.psew, bl.comew, bl.tg, bl.income, bl.work_lastmonth, bl.mue, bl.rew, bl.vew,
                      bl.rm, bl.rmasp, bl.rmtime, bl.ss, blp.ss, blf.ss, bl.lead, bl.goal, bl.log.exp.pc) %>%
  group_by(workshop_name) %>% summarize(sd.age = sd(age, na.rm=T),
                                        sd.hhsize = sd(hhsize, na.rm=T), 
                                        sd.sa = sd(bl.saew, na.rm=T),
                                        sd.ps = sd(bl.psew, na.rm=T),
                                        sd.nlead = sd(bl.nlead, na.rm=T),
                                        sd.conhv = sd(bl.conhvew, na.rm=T),
                                        sd.con = sd(bl.conew, na.rm=T),
                                        sd.com = sd(bl.comew, na.rm=T),
                                        sd.lead.ind = sd(bl.lead.ind, na.rm=T),
                                        sd.ss = sd(bl.ss, na.rm=T),
                                        sd.pss = sd(blp.ss, na.rm=T),
                                        sd.fss = sd(blf.ss, na.rm=T),
                                        sd.v = sd(bl.vew, na.rm=T),
                                        sd.r = sd(bl.rew, na.rm=T),
                                        sd.lead = sd(bl.lead, na.rm=T),
                                        sd.tg = sd(bl.tg, na.rm=T),
                                        sd.rm = sd(bl.rm, na.rm=T),
                                        sd.rmasp = sd(bl.rmasp, na.rm=T),
                                        sd.rmtime = sd(bl.rmtime, na.rm=T),
                                        sd.goal = sd(bl.goal, na.rm=T),
                                        sd.income = sd(bl.income, na.rm=T), 
                                        sd.work = sd(bl.work_lastmonth, na.rm=T),
                                        sd.mue = sd(bl.mue, na.rm=T),
                                        sd.cons = sd(bl.log.exp.pc, na.rm=T))

balance.means <- as.data.frame(t(as.matrix(balance.means)))
colnames(balance.means) <- balance.means[1,]
balance.means <- balance.means[c(2:25),]
balance.means$Control <- as.numeric(balance.means$Control)
balance.means$`Professional Development` <- as.numeric(balance.means$`Professional Development`)
balance.means$`Storytelling for Leadership` <- as.numeric(balance.means$`Storytelling for Leadership`)
balance.means <- balance.means %>% mutate_if(is.numeric, round, digits=2)

balance.sds <- as.data.frame(t(as.matrix(balance.sds)))
colnames(balance.sds) <- balance.sds[1,]
balance.sds <- balance.sds[c(2:25),]
balance.sds$Control <- as.numeric(balance.sds$Control)
balance.sds$`Professional Development` <- as.numeric(balance.sds$`Professional Development`)
balance.sds$`Storytelling for Leadership` <- as.numeric(balance.sds$`Storytelling for Leadership`)
balance.sds <- balance.sds %>% mutate_if(is.numeric, round, digits=2)

table <- matrix(NA, nrow=28, ncol=5)

table[1,] <- c("","Cash Control", "Soft Skills", "Storytelling","p-value")
table[2,] <- unlist(c("Age", balance.means[1,1:3],round(unlist(summary(bal.age, robust = TRUE)[11])[4],2)))
table[3,] <- c("",paste0("(",balance.sds[1,1],")"),paste0("(",balance.sds[1,2],")"),paste0("(",balance.sds[1,3],")"),"")
table[4,] <- unlist(c("HH Size", balance.means[2,1:3],round(unlist(summary(bal.hhsize, robust = TRUE)[11])[4],2)))
table[5,] <- c("",paste0("(",balance.sds[2,1],")"),paste0("(",balance.sds[2,2],")"),paste0("(",balance.sds[2,3],")"),"")
table[6,] <- unlist(c("Self-Advocacy", balance.means[3,1:3],round(unlist(summary(bal.sa, robust = TRUE)[11])[4],2)))
table[7,] <- c("",paste0("(",balance.sds[3,1],")"),paste0("(",balance.sds[3,2],")"),paste0("(",balance.sds[3,3],")"),"")
table[8,] <- unlist(c("Problem Solving", balance.means[4,1:3],round(unlist(summary(bal.ps, robust = TRUE)[11])[4],2)))
table[9,] <- c("",paste0("(",balance.sds[4,1],")"),paste0("(",balance.sds[4,2],")"),paste0("(",balance.sds[4,3],")"),"")
table[10,] <- unlist(c("No. Leadership Positions", balance.means[5,1:3],round(unlist(summary(bal.nlead, robust = TRUE)[11])[4],2)))
table[11,] <- c("",paste0("(",balance.sds[5,1],")"),paste0("(",balance.sds[5,2],")"),paste0("(",balance.sds[5,3],")"),"")
table[12,] <- unlist(c("Connectedness to Home/Village", balance.means[6,1:3],round(unlist(summary(bal.conhv, robust = TRUE)[11])[4],2)))
table[13,] <- c("",paste0("(",balance.sds[6,1],")"),paste0("(",balance.sds[6,2],")"),paste0("(",balance.sds[6,3],")"),"")
table[14,] <- unlist(c("Connectedness to Co-op", balance.means[7,1:3],round(unlist(summary(bal.con, robust = TRUE)[11])[4],2)))
table[15,] <- c("",paste0("(",balance.sds[7,1],")"),paste0("(",balance.sds[7,2],")"),paste0("(",balance.sds[7,3],")"),"")
table[16,] <- unlist(c("Communication", balance.means[8,1:3],round(unlist(summary(bal.com, robust = TRUE)[11])[4],2)))
table[17,] <- c("",paste0("(",balance.sds[8,1],")"),paste0("(",balance.sds[8,2],")"),paste0("(",balance.sds[8,3],")"),"")
table[18,] <- unlist(c("Any Leadership Position", balance.means[9,1:3],round(unlist(summary(bal.lead.ind, robust = TRUE)[11])[4],2)))
table[19,] <- c("",paste0("(",balance.sds[9,1],")"),paste0("(",balance.sds[9,2],")"),paste0("(",balance.sds[9,3],")"),"")
table[20,] <- unlist(c("Cash Earnings", balance.means[21,1:3],round(unlist(summary(bal.inc, robust = TRUE)[11])[4],2)))
table[21,] <- c("",paste0("(",balance.sds[21,1],")"),paste0("(",balance.sds[21,2],")"),paste0("(",balance.sds[21,3],")"),"")
table[22,] <- unlist(c("Earn Any Cash", balance.means[22,1:3],round(unlist(summary(bal.inc, robust = TRUE)[11])[4],2)))
table[23,] <- c("",paste0("(",balance.sds[22,1],")"),paste0("(",balance.sds[22,2],")"),paste0("(",balance.sds[22,3],")"),"")
table[24,] <- unlist(c("MUE", balance.means[23,1:3],round(unlist(summary(bal.mue, robust = TRUE)[11])[4],2)))
table[25,] <- c("",paste0("(",balance.sds[23,1],")"),paste0("(",balance.sds[23,2],")"),paste0("(",balance.sds[23,3],")"),"")
table[26,] <- unlist(c("Log(Expenditures per Capita)", balance.means[24,1:3],round(unlist(summary(bal.cons, robust = TRUE)[11])[4],2)))
table[27,] <- c("",paste0("(",balance.sds[24,1],")"),paste0("(",balance.sds[24,2],")"),paste0("(",balance.sds[24,3],")"),"")
table[28,] <- unlist(c("N",length(bl$resp_id[which(bl$workshop_name=="Control")]),length(bl$resp_id[which(bl$workshop_name=="Professional Development")]),
                       length(bl$resp_id[which(bl$workshop_name=="Storytelling for Leadership")]),""))

myxtable <- xtable(table)

align(myxtable) <- c("l",rep("c",5))

caption(myxtable) <- "Balance"

label(myxtable) <- "balance"

print(myxtable, floating = TRUE, caption.placement="top",sanitize.text.function = identity,
      include.colnames=F,include.rownames=F,table.placement="H",
      hline.after=NULL,
      add.to.row=list(
        pos = list(0,1,nrow(table)-1,nrow(table)),
        command = c(
          paste0("\\toprule"),
          " \\midrule ",
          " \\midrule ",
          "\\bottomrule \\multicolumn{5}{l}{\\multirow{2}{14cm}{\\small Notes: Mean baseline covariates by treatment group. Standard deviations are in parentheses. Column 4 reports p-values associated with F-tests of joint equality between the three groups.}} \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ "
        )),
      type = "latex",file="balance.tex")

##########################################################################################
#Appendix Figures and Tables  
##########################################################################################
#####Figure A1 - Distribution of Women from Co-ops#####
#Number of women in our sample from the same cooperative
c.count <- mutate(full, counter = 1) %>% group_by(Coop) %>% summarize(n = sum(counter))

median(c.count$n)

coop.dist <- ggplot(data = c.count, aes(n)) + geom_histogram(binwidth = 1) + 
  labs(x = "Number of Women", 
       y = "Number of Cooperatives") + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5))

ggsave(filename="../Figures/Coops.eps", plot=coop.dist, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)

#####Figure A2 - Distribution of baseline monthly cash earnings####
inc.dist <- ggplot(data = full[which(full$bl.income<150000),], aes(bl.income)) + 
  geom_histogram(binwidth = 5000) + 
  geom_vline(xintercept = median(full$bl.income, na.rm=T), color = "grey", linetype="dashed") + 
  geom_vline(xintercept = mean(full$bl.income, na.rm=T), color = "grey") + 
  geom_vline(xintercept = 30000, color = "black") + 
  labs(x = "Income earned the previous month at baseline", 
       y = "No. Participants") + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), 
                     axis.text.x = element_text(angle = 45, hjust = 1)) 

ggsave(filename="../Figures/Inc_Dist.eps", plot=inc.dist, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)

#####Figure A3 - Treatment effect - speaking up#####
su.sd <- sqrt(mean(full$fu.su[which(full$workshop_name=="Control")],na.rm=T)*(1-mean(full$fu.su[which(full$workshop_name=="Control")],na.rm=T)))
su.pd <- as.data.frame(cbind(su.reg$coefficients, su.reg$rse))
su.pd <- cbind(su.pd, rownames(su.pd))
colnames(su.pd) <- c("ATE","RSE","Tmt")
su.pd <- mutate(su.pd, dep.var = "Speaking Up")

su.pd <- filter(su.pd, Tmt=="SFL"|Tmt=="PD") %>%
  mutate(SD = su.sd,
         ATE.std = ATE/SD, RSE.std = RSE/SD, CI.high = ATE.std + 1.96*RSE.std,
         CI.low = ATE.std - 1.96*RSE.std)
su.pd <- mutate(su.pd, Treatment = if_else(Tmt=="PD","Soft Skills","Storytelling"))

su.plot <- ggplot(data = su.pd, aes(y = ATE.std, x = Treatment, group = Treatment,
                                    fill = Treatment, color = Treatment, 
                                    ymin = CI.low, ymax = CI.high)) + 
  geom_point(position = position_dodge(0.9)) + 
  geom_errorbar(width = 0.1, position = position_dodge(0.9)) + 
  geom_hline(yintercept = 0, color = "black") + 
  labs(x = "Treatment", 
       y = "Treatment Effect \n (Standard deviations)") + 
  scale_fill_grey(start = 0.5, end = 0.05) + scale_color_grey(start = 0.5, end = 0.05) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), 
                     axis.text.x = element_text(angle = 45, hjust = 1))

ggsave(filename="../Figures/SU.eps", plot=su.plot, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)
#####Table A1 - Primary Psychosocial ATEs in levels ####
SFL.pps$adj.SE <- round(SFL.pps$adj.SE,2)
PD.pps$adj.SE <- round(PD.pps$adj.SE,2)

table <- matrix(NA, nrow=13, ncol=9)

table[1,] <- c("","Self-Advocacy","Problem Solving","Leadership Positions","Connection Home/Village","Connection Co-op",
               "Communication Skills","Any Leadership","Index")
table[2,] <- unlist(c("Storytelling", round(sa.reg.ew$coefficients[2],2), round(ps.reg.ew$coefficients[2],2), round(nlead.reg$coefficients[2],2), 
                      round(conhv.reg.ew$coefficients[2],2), round(con.reg.ew$coefficients[2],2), round(com.reg.ew$coefficients[2],2), 
                      round(lead.ind.reg$coefficients[2],2), round(pps.reg$coefficients[2],2)))
table[3,] <- c("",paste0("(",round(sa.reg.ew$rse[2],2),")"), paste0("(",round(ps.reg.ew$rse[2],2),")"), paste0("(",round(nlead.reg$rse[2],2),")"),
               paste0("(",round(conhv.reg.ew$rse[2],2),")"), paste0("(",round(con.reg.ew$rse[2],2),")"), paste0("(",round(com.reg.ew$rse[2],2),")"),
               paste0("(",round(lead.ind.reg$rse[2],2),")"), paste0("(",round(pps.reg$rse[2],2),")"))
table[4,] <- c("",paste0("(",SFL.pps$adj.SE[3],")"), paste0("(",SFL.pps$adj.SE[1],")"), paste0("(",SFL.pps$adj.SE[6],")"),
               paste0("(",SFL.pps$adj.SE[2],")"), paste0("(",SFL.pps$adj.SE[7],")"), paste0("(",SFL.pps$adj.SE[4],")"),
               paste0("(",SFL.pps$adj.SE[5],")"), "")
table[5,] <- unlist(c("Soft Skills", round(sa.reg.ew$coefficients[3],2), round(ps.reg.ew$coefficients[3],2), round(nlead.reg$coefficients[3],2), 
                      round(conhv.reg.ew$coefficients[3],2), round(con.reg.ew$coefficients[3],2), round(com.reg.ew$coefficients[3],2), 
                      round(lead.ind.reg$coefficients[3],2), round(pps.reg$coefficients[3],2)))
table[6,] <- c("",paste0("(",round(sa.reg.ew$rse[3],2),")"), paste0("(",round(ps.reg.ew$rse[3],2),")"), paste0("(",round(nlead.reg$rse[3],2),")"),
               paste0("(",round(conhv.reg.ew$rse[3],2),")"), paste0("(",round(con.reg.ew$rse[3],2),")"), paste0("(",round(com.reg.ew$rse[3],2),")"),
               paste0("(",round(lead.ind.reg$rse[3],2),")"), paste0("(",round(pps.reg$rse[3],2),")"))
table[7,] <- c("",paste0("(",PD.pps$adj.SE[3],")"), paste0("(",PD.pps$adj.SE[4],")"), paste0("(",PD.pps$adj.SE[6],")"),
               paste0("(",PD.pps$adj.SE[2],")"), paste0("(",PD.pps$adj.SE[1],")"), paste0("(",PD.pps$adj.SE[7],")"),
               paste0("(",PD.pps$adj.SE[5],")"), "")
table[8,] <- unlist(c("Baseline", round(sa.reg.ew$coefficients[4],2), round(ps.reg.ew$coefficients[4],2), round(nlead.reg$coefficients[4],2),
                      round(conhv.reg.ew$coefficients[4],2), round(con.reg.ew$coefficients[4],2), round(com.reg.ew$coefficients[4],2), 
                      round(lead.ind.reg$coefficients[4],2), round(pps.reg$coefficients[4],2)))
table[9,] <- c("Level",paste0("(",round(sa.reg.ew$rse[4],2),")"), paste0("(",round(ps.reg.ew$rse[4],2),")"), paste0("(",round(nlead.reg$rse[4],2),")"),
               paste0("(",round(conhv.reg.ew$rse[4],2),")"), paste0("(",round(con.reg.ew$rse[4],2),")"), paste0("(",round(com.reg.ew$rse[4],2),")"),
               paste0("(",round(lead.ind.reg$rse[4],2),")"), paste0("(",round(pps.reg$rse[4],2),")"))
table[10,] <- unlist(c("Control Mean",round(sa.mean,2), round(ps.mean,2), round(l.mean,2), round(chv.mean,2), round(cc.mean,2), 
                       round(com.mean,2), round(lind.mean,2), round(pps.mean,2)))
table[11,] <- unlist(c("N",sa.reg.ew$N, ps.reg.ew$N, nlead.reg$N, conhv.reg.ew$N, con.reg.ew$N, com.reg.ew$N, lead.ind.reg$N, pps.reg$N))
table[12,] <- c("Adj. R2",round(as.numeric(summary(sa.reg.ew)[15]),2),round(as.numeric(summary(ps.reg.ew)[15]),2),
                round(as.numeric(summary(nlead.reg)[15]),2), round(as.numeric(summary(conhv.reg.ew)[15]),2), 
                round(as.numeric(summary(con.reg.ew)[15]),2),round(as.numeric(summary(com.reg.ew)[15]),2),
                round(as.numeric(summary(lead.ind.reg)[15]),2),round(as.numeric(summary(pps.reg)[15]),2))
table[13,] <- unlist(c("S = SS (p-value)", round(linearHypothesis(sa.reg.ew, "SFL = PD")[2,4],2), round(linearHypothesis(ps.reg.ew, "SFL = PD")[2,4],2),
                       round(linearHypothesis(nlead.reg, "SFL = PD")[2,4],2), round(linearHypothesis(conhv.reg.ew, "SFL = PD")[2,4],2), round(linearHypothesis(con.reg.ew, "SFL = PD")[2,4],2),
                       round(linearHypothesis(com.reg.ew, "SFL = PD")[2,4],2), round(linearHypothesis(lead.ind.reg, "SFL = PD")[2,4],2), round(linearHypothesis(pps.reg, "SFL = PD")[2,4],2)))

myxtable <- xtable(table)

align(myxtable) <- c("l","p{2cm}", rep("p{2.1cm}",3), rep("p{2.4cm}",4),"p{2cm}")

caption(myxtable) <- "Primary Psychosocial Outcomes"

label(myxtable) <- "ppstab"

print(myxtable, floating = TRUE, caption.placement="top",sanitize.text.function = identity,
      include.colnames=F,include.rownames=F,table.placement="H",
      hline.after=NULL,
      add.to.row=list(
        pos = list(0,1,9,nrow(table)),
        command = c(
          paste0("\\toprule"),
          " \\midrule ",
          " \\midrule ",
          "\\bottomrule \\multicolumn{9}{l}{\\multirow{2}{22.5cm}{\\small Notes: Average treatment effects for our family of pre-registered, primary psychosocial outcomes. Self-advocacy, problem-solving, connection to home/village, connection to co-op, and communication skills are pre-registered indices, typically with all questions on 5-point Likert scales. Leadership positions is the number of leadership positions currently held and leadership position is an indicator variable equal to one if the respondent currently holds any leadership positions.  Index is an index of all seven outcomes (note that we did not pre-register this index). We show naive White robust standard errors in the first set of parentheses. The second set of parentheses shows standard errors adjusted to control the false discovery rate within this family of outcomes according to the step-up procedure in Benjamini, Kreiger, and Yekutieli (2006). The bottom row shows the p-values associated with hypothesis tests of equality between the estimated average treatment effect for Storytelling (S) versus Soft Skills (SS).}} \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ "
        )),
      type = "latex",file="ppstable.tex")

#####Table A2 - Secondary Psychosocial ATEs in levels ####
SFL.sps$adj.SE <- round(SFL.sps$adj.SE,2)
PD.sps$adj.SE <- round(PD.sps$adj.SE,2)

table <- matrix(NA, nrow=13, ncol=14)

table[1,] <- c("","SSS","Past SSS","Future SSS", "CWG","Self-value","Resilient","ID Lead","Role Model",
               "Aspire RM","Time RM","Has Goal","Achieve Goal","Index")
table[2,] <- unlist(c("Story-", round(ss.reg$coefficients[2],2), round(ssp.reg$coefficients[2],2), round(ssf.reg$coefficients[2],2), 
                      round(tg.reg$coefficients[2],2), round(val.reg.ew$coefficients[2],2), round(res.reg.ew$coefficients[2],2), 
                      round(lead.reg$coefficients[2],2), round(rm.reg$coefficients[2],2), round(rmasp.reg$coefficients[2],2),
                      round(rmtime.reg$coefficients[2],2), round(goal.reg$coefficients[2],2), round(com.goal.reg$coefficients[2],2),
                      round(sps.reg$coefficients[2],2)))
table[3,] <- c("telling",paste0("(",round(ss.reg$rse[2],2),")"), paste0("(",round(ssp.reg$rse[2],2),")"), paste0("(",round(ssf.reg$rse[2],2),")"),
               paste0("(",round(tg.reg$rse[2],2),")"), paste0("(",round(val.reg.ew$rse[2],2),")"), paste0("(",round(res.reg.ew$rse[2],2),")"),
               paste0("(",round(lead.reg$rse[2],2),")"), paste0("(",round(rm.reg$rse[2],2),")"), paste0("(",round(rmasp.reg$rse[2],2),")"),
               paste0("(",round(rmtime.reg$rse[2],2),")"), paste0("(",round(goal.reg$rse[2],2),")"), paste0("(",round(com.goal.reg$rse[2],2),")"),
               paste0("(",round(sps.reg$rse[2],2),")"))
table[4,] <- c("",paste0("(",SFL.sps$adj.SE[3],")"), paste0("(",SFL.sps$adj.SE[2],")"), paste0("(",SFL.sps$adj.SE[10],")"),
               paste0("(",SFL.sps$adj.SE[12],")"), paste0("(",SFL.sps$adj.SE[4],")"), paste0("(",SFL.sps$adj.SE[8],")"),
               paste0("(",SFL.sps$adj.SE[1],")"),paste0("(",SFL.sps$adj.SE[5],")"),paste0("(",SFL.sps$adj.SE[7],")"),
               paste0("(",SFL.sps$adj.SE[9],")"),paste0("(",SFL.sps$adj.SE[11],")"),paste0("(",SFL.sps$adj.SE[6],")"),"")
table[5,] <- unlist(c("Soft", round(ss.reg$coefficients[3],2), round(ssp.reg$coefficients[3],2), round(ssf.reg$coefficients[3],2), 
                      round(tg.reg$coefficients[3],2), round(val.reg.ew$coefficients[3],2), round(res.reg.ew$coefficients[3],2), 
                      round(lead.reg$coefficients[3],2), round(rm.reg$coefficients[3],2), round(rmasp.reg$coefficients[3],2),
                      round(rmtime.reg$coefficients[3],2), round(goal.reg$coefficients[3],2), round(com.goal.reg$coefficients[3],2),
                      round(sps.reg$coefficients[3],2)))
table[6,] <- c("Skills",paste0("(",round(ss.reg$rse[3],2),")"), paste0("(",round(ssp.reg$rse[3],2),")"), paste0("(",round(ssf.reg$rse[3],2),")"),
               paste0("(",round(tg.reg$rse[3],2),")"), paste0("(",round(val.reg.ew$rse[3],2),")"), paste0("(",round(res.reg.ew$rse[3],2),")"),
               paste0("(",round(lead.reg$rse[3],2),")"), paste0("(",round(rm.reg$rse[3],2),")"), paste0("(",round(rmasp.reg$rse[3],2),")"),
               paste0("(",round(rmtime.reg$rse[3],2),")"), paste0("(",round(goal.reg$rse[3],2),")"), paste0("(",round(com.goal.reg$rse[3],2),")"),
               paste0("(",round(sps.reg$rse[3],2),")"))
table[7,] <- c("",paste0("(",PD.sps$adj.SE[11],")"), paste0("(",PD.sps$adj.SE[6],")"), paste0("(",PD.sps$adj.SE[4],")"),
               paste0("(",PD.sps$adj.SE[5],")"), paste0("(",PD.sps$adj.SE[3],")"), paste0("(",PD.sps$adj.SE[12],")"),
               paste0("(",PD.sps$adj.SE[1],")"),paste0("(",PD.sps$adj.SE[9],")"),paste0("(",PD.sps$adj.SE[7],")"),
               paste0("(",PD.sps$adj.SE[10],")"),paste0("(",PD.sps$adj.SE[8],")"),paste0("(",PD.sps$adj.SE[2],")"),"")
table[8,] <- unlist(c("Baseline", round(ss.reg$coefficients[4],2), round(ssp.reg$coefficients[4],2), round(ssf.reg$coefficients[4],2), 
                      round(tg.reg$coefficients[4],2), round(val.reg.ew$coefficients[4],2), round(res.reg.ew$coefficients[4],2), 
                      round(lead.reg$coefficients[4],2), round(rm.reg$coefficients[4],2), round(rmasp.reg$coefficients[4],2),
                      round(rmtime.reg$coefficients[4],2), round(goal.reg$coefficients[4],2), round(com.goal.reg$coefficients[4],2),
                      round(sps.reg$coefficients[4],2)))
table[9,] <- c("Level",paste0("(",round(ss.reg$rse[4],2),")"), paste0("(",round(ssp.reg$rse[4],2),")"), paste0("(",round(ssf.reg$rse[4],2),")"),
               paste0("(",round(tg.reg$rse[4],2),")"), paste0("(",round(val.reg.ew$rse[4],2),")"), paste0("(",round(res.reg.ew$rse[4],2),")"),
               paste0("(",round(lead.reg$rse[4],2),")"), paste0("(",round(rm.reg$rse[4],2),")"), paste0("(",round(rmasp.reg$rse[4],2),")"),
               paste0("(",round(rmtime.reg$rse[4],2),")"), paste0("(",round(goal.reg$rse[4],2),")"), paste0("(",round(com.goal.reg$rse[4],2),")"),
               paste0("(",round(sps.reg$rse[4],2),")"))
table[10,] <- unlist(c("Control Mean",round(ss.mean,2), round(pss.mean,2), round(fss.mean,2), round(conw.mean,2), round(v.mean,2), 
                       round(r.mean,2), round(idl.mean,2), round(rm.mean,2), round(rma.mean,2), round(rmt.mean,2), round(g.mean,2), 
                       round(gc.mean,2), round(sps.mean,2)))
table[11,] <- unlist(c("N",ss.reg$N, ssp.reg$N, ssf.reg$N, tg.reg$N, val.reg.ew$N, res.reg.ew$N, lead.reg$N, rm.reg$N, rmasp.reg$N, rmtime.reg$N,
                       goal.reg$N, com.goal.reg$N, sps.reg$N))
table[12,] <- unlist(c("Adj. R2",round(as.numeric(summary(ss.reg)[15]),2),round(as.numeric(summary(ssp.reg)[15]),2),
                       round(as.numeric(summary(ssf.reg)[15]),2), round(as.numeric(summary(tg.reg)[15]),2), 
                       round(as.numeric(summary(val.reg.ew)[15]),2),round(as.numeric(summary(res.reg.ew)[15]),2),
                       round(as.numeric(summary(lead.reg)[15]),2),round(as.numeric(summary(rm.reg)[15]),2),
                       round(as.numeric(summary(rmasp.reg)[15]),2), round(as.numeric(summary(rmtime.reg)[15]),2),
                       round(as.numeric(summary(goal.reg)[15]),2),round(as.numeric(summary(com.goal.reg)[15]),2),
                       round(as.numeric(summary(sps.reg)[15]),2)))
table[13,] <- unlist(c("S = SS", round(linearHypothesis(ss.reg, "SFL = PD")[2,4],2), round(linearHypothesis(ssp.reg, "SFL = PD")[2,4],2),
                       round(linearHypothesis(ssf.reg, "SFL = PD")[2,4],2), round(linearHypothesis(tg.reg, "SFL = PD")[2,4],2), 
                       round(linearHypothesis(val.reg.ew, "SFL = PD")[2,4],2), round(linearHypothesis(res.reg.ew, "SFL = PD")[2,4],2), 
                       round(linearHypothesis(lead.reg, "SFL = PD")[2,4],2), round(linearHypothesis(rm.reg, "SFL = PD")[2,4],2),
                       round(linearHypothesis(rmasp.reg, "SFL = PD")[2,4],2),round(linearHypothesis(rmtime.reg, "SFL = PD")[2,4],2),
                       round(linearHypothesis(goal.reg, "SFL = PD")[2,4],2),round(linearHypothesis(com.goal.reg, "SFL = PD")[2,4],2),
                       round(linearHypothesis(sps.reg, "SFL = PD")[2,4],2)))

myxtable <- xtable(table)

align(myxtable) <- c("l","p{1.3cm}",rep("p{1.2cm}",13))

caption(myxtable) <- "Secondary Psychosocial Outcomes"

label(myxtable) <- "spstab"

print(myxtable, floating = TRUE, caption.placement="top",sanitize.text.function = identity,
      include.colnames=F,include.rownames=F,table.placement="H",
      hline.after=NULL,
      add.to.row=list(
        pos = list(0,1,9,nrow(table)),
        command = c(
          paste0("\\toprule"),
          " \\midrule ",
          " \\midrule ",
          "\\bottomrule \\multicolumn{14}{l}{\\multirow{2}{22.5cm}{\\small Notes: Average treatment effects for our family of pre-registered, secondary psychosocial outcomes. SSS is current subjective social status. SSS Past is SSS evaluated 5 years in the past and SSS Future is expected SSS five years in the future. CWG (connection to workshop group), self-value, and resilience are pre-registered indices of questions. ID lead measures how closely the respondent identifies with being a leader. Role model is an indicator variable equal to one if the respondent has a role model. Aspire RM is an indicator equal to one if the respondent aspires to be their role model. Time RM is the time (in years) that the respondent thinks it will take them to be like their role model. Goal is an indicator equal to one if the respondent has a goal at follow-up. Achieve goal is an indicator equal to one if the respondent reports that they achieved the goal set at baseline by the time of the follow-up survey. Index is an index of all outcomes except aspiring to be the role model and the time to take to be like the role model, to avoid losses in power (note that we did not pre-register this index). We show naive White robust standard errors in the first set of parentheses. The second set of parentheses shows standard errors adjust to control the false discovery rate within this family of outcomes according to the step-up procedure in Benjamini, Kreiger, and Yekutieli (2006). The bottom row shows the p-values associated with hypothesis tests of equality between the estimated average treatment effect for Storytelling (S) versys Soft Skills (SS).}} \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ "
        )),
      type = "latex",file="spstable.tex")
#####Table A3 - Economic ATEs in levels and non-transformed ####
#Estimate effects for consumption spending in levels
cons <- felm(fu.total.exp ~ SFL + PD + bl.total.exp |0|0, data = full)
summary(cons, robust = TRUE)
linearHypothesis(cons, "SFL = PD")

SFL.e$adj.SE <- round(SFL.e$adj.SE,2)
PD.e$adj.SE <- round(PD.e$adj.SE,2)

table <- matrix(NA, nrow=13, ncol=7)

table[1,] <- c("","Marginal Utility of Expenditure","Log(Consumption Expenditures per Capita)","Consumption Expenditures per Capita (RWF)","Monthly Cash Earnings (RWF)","Earn Any Cash","Index")
table[2,] <- unlist(c("Storytelling", round(mue.reg$coefficients[2],2), round(exp.reg$coefficients[2],2), round(cons$coefficients[2],2),
                      round(inc$coefficients[2],2), round(work.reg$coefficients[2],2),
                      round(e.reg$coefficients[2],2)))
table[3,] <- c("",paste0("(",round(mue.reg$rse[2],2),")"), paste0("(",round(exp.reg$rse[2],2),")"),paste0("(",round(cons$rse[2],2),")"),
               paste0("(",round(inc$rse[2],2),")"), paste0("(",round(work.reg$rse[2],2),")"),
               paste0("(",round(e.reg$rse[2],2),")"))
table[4,] <- c("",paste0("(",SFL.e$adj.SE[3],")"), paste0("(",SFL.e$adj.SE[4],")"), "", paste0("(",SFL.e$adj.SE[2],")"), paste0("(",SFL.e$adj.SE[1],")"),"")
table[5,] <- unlist(c("Soft Skills", round(mue.reg$coefficients[3],2), round(exp.reg$coefficients[3],2), round(cons$coefficients[3],2),
                      round(inc$coefficients[3],2), round(work.reg$coefficients[3],2),
                      round(e.reg$coefficients[3],2)))
table[6,] <- c("",paste0("(",round(mue.reg$rse[3],2),")"), paste0("(",round(exp.reg$rse[3],2),")"), paste0("(",round(cons$rse[3],2),")"),
               paste0("(",round(inc$rse[3],2),")"), paste0("(",round(work.reg$rse[3],2),")"),
               paste0("(",round(e.reg$rse[3],2),")"))
table[7,] <- c("",paste0("(",PD.e$adj.SE[3],")"), paste0("(",PD.e$adj.SE[4],")"),"", paste0("(",PD.e$adj.SE[2],")"), paste0("(",PD.e$adj.SE[1],")"),"")
table[8,] <- unlist(c("Baseline Level", round(mue.reg$coefficients[4],2), round(cons$coefficients[4],2),round(exp.reg$coefficients[4],2),
                      round(inc$coefficients[4],2), round(work.reg$coefficients[4],2),
                      round(e.reg$coefficients[4],2)))
table[9,] <- c("",paste0("(",round(mue.reg$rse[4],2),")"), paste0("(",round(exp.reg$rse[4],2),")"), paste0("(",round(cons$rse[4],2),")"),
               paste0("(",round(inc$rse[4],2),")"), paste0("(",round(work.reg$rse[4],2),")"),
               paste0("(",round(e.reg$rse[4],2),")"))
table[10,] <- unlist(c("Control Mean",round(mue.mean,2), round(exp.mean,2), round(mean(full$fu.exp.pc[which(full$workshop_name=="Control")],na.rm=T),2), round(inc.mean,2), round(inci.mean,2),round(e.mean,2)))
table[11,] <- unlist(c("N",mue.reg$N, exp.reg$N, cons$N , inc$N, work.reg$N, e.reg$N))
table[12,] <- c("Adj. R2",round(as.numeric(summary(mue.reg)[15]),2), round(as.numeric(summary(exp.reg)[15]),2),round(as.numeric(summary(cons)[15]),2), round(as.numeric(summary(inc)[15]),2),
                round(as.numeric(summary(work.reg)[15]),2), round(as.numeric(summary(e.reg)[15]),2))
table[13,] <- unlist(c("S = SS (p-value)", round(linearHypothesis(mue.reg, "SFL = PD")[2,4],2), round(linearHypothesis(exp.reg, "SFL = PD")[2,4],2),round(linearHypothesis(cons, "SFL = PD")[2,4],2),
                       round(linearHypothesis(inc, "SFL = PD")[2,4],2),
                       round(linearHypothesis(work.reg, "SFL = PD")[2,4],2), round(linearHypothesis(e.reg, "SFL = PD")[2,4],2)))

myxtable <- xtable(table)

align(myxtable) <- c("l","p{4cm}",rep("p{2.7cm}",6))

caption(myxtable) <- "Economic Outcomes"

label(myxtable) <- "etab"

print(myxtable, floating = TRUE, caption.placement="top",sanitize.text.function = identity,
      include.colnames=F,include.rownames=F,table.placement="H",
      hline.after=NULL,
      add.to.row=list(
        pos = list(0,1,9,nrow(table)),
        command = c(
          paste0("\\toprule"),
          " \\midrule ",
          " \\midrule ",
          "\\bottomrule \\multicolumn{7}{l}{\\multirow{2}{22.5cm}{\\small Notes: Average treatment effects for our family of pre-registered economic outcomes, plus per capita consumption expenditures on 18 food items in logs and levels.  Index is an index of the marginal utility of expenditure, monthly cash earnings, and the log of per capita consumption expenditure on the 18 food items that we use to compute the marginal utility of expenditure (note that we did not pre-register this index). We show naive White robust standard errors in the first set of parentheses. The second set of parentheses shows standard errors adjust to control the false discovery rate within this family of outcomes according to the step-up procedure in Benjamini, Kreiger, and Yekutieli (2006). The bottom row shows the p-values associated with hypothesis tests of equality between the estimated average treatment effect for Storytelling (S) versus Soft Skills (SS).}} \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ "
        )),
      type = "latex",file="etable.tex")


#####Table A4 - Spillovers#####
table <- matrix(NA, nrow=14, ncol=7)

table[1,] <- c("","(1)","(2)","(3)","(4)","(5)","(6)")
table[2,] <- unlist(c("No. Storytelling", round(c.pps.reg$coefficients[2],2),"", round(c.sps.reg$coefficients[2],2), "",
                      round(c.e.reg$coefficients[2],2),""))
table[3,] <- c("",paste0("(",round(c.pps.reg$rse[2],2),")"), "", paste0("(",round(c.sps.reg$rse[2],2),")"), "",
               paste0("(",round(c.e.reg$rse[2],2),")"),"")
table[4,] <- unlist(c("No. Soft Skills", round(c.pps.reg$coefficients[3],2), "", round(c.sps.reg$coefficients[3],2), "",
                      round(c.e.reg$coefficients[3],2),""))
table[5,] <- c("",paste0("(",round(c.pps.reg$rse[3],2),")"), "", paste0("(",round(c.sps.reg$rse[3],2),")"), "",
               paste0("(",round(c.e.reg$rse[3],2),")"),"")
table[6,] <- unlist(c("Above-Median Storytelling", "",round(c2.pps.reg$coefficients[2],2), "", round(c2.sps.reg$coefficients[2],2), "",
                      round(c2.e.reg$coefficients[2],2)))
table[7,] <- c("", "", paste0("(",round(c2.pps.reg$rse[2],2),")"), "", paste0("(",round(c2.sps.reg$rse[2],2),")"), "",
               paste0("(",round(c2.e.reg$rse[2],2),")"))
table[8,] <- unlist(c("Above-Median Soft Skills", "", round(c2.pps.reg$coefficients[3],2), "", round(c2.sps.reg$coefficients[3],2), "",
                      round(c2.e.reg$coefficients[3],2)))
table[9,] <- c("","",paste0("(",round(c2.pps.reg$rse[3],2),")"), "", paste0("(",round(c2.sps.reg$rse[3],2),")"), "",
               paste0("(",round(c2.e.reg$rse[3],2),")"))
table[10,] <- unlist(c("Baseline Level", round(c.pps.reg$coefficients[1],2), round(c.sps.reg$coefficients[1],2),
                       round(c.e.reg$coefficients[1],2), round(c2.pps.reg$coefficients[1],2), round(c2.sps.reg$coefficients[1],2),
                       round(c2.e.reg$coefficients[1],2)))
table[11,] <- c("",paste0("(",round(c.pps.reg$rse[1],2),")"), paste0("(",round(c.sps.reg$rse[1],2),")"),
                paste0("(",round(c.e.reg$rse[1],2),")"), paste0("(",round(c2.pps.reg$rse[1],2),")"), paste0("(",round(c2.sps.reg$rse[1],2),")"),
                paste0("(",round(c2.e.reg$rse[1],2),")"))

table[12,] <- unlist(c("N",c.pps.reg$N, c.sps.reg$N, c.e.reg$N, c2.pps.reg$N, c2.sps.reg$N, c2.e.reg$N))
table[13,] <- c("Adj. R2",round(as.numeric(summary(c.pps.reg)[15]),2), round(as.numeric(summary(c.sps.reg)[15]),2),
                round(as.numeric(summary(c.e.reg)[15]),2),round(as.numeric(summary(c2.pps.reg)[15]),2), round(as.numeric(summary(c2.sps.reg)[15]),2),
                round(as.numeric(summary(c2.e.reg)[15]),2))
table[14,] <- unlist(c("S=SS (p-value)", round(linearHypothesis(c.pps.reg, "c.SFL = c.pd")[2,4],2), 
                       round(linearHypothesis(c.sps.reg, "c.SFL = c.pd")[2,4],2),
                       round(linearHypothesis(c.e.reg, "c.SFL = c.pd")[2,4],2),round(linearHypothesis(c2.pps.reg, "c.sfl.high = c.pd.high")[2,4],2), 
                       round(linearHypothesis(c2.sps.reg, "c.sfl.high = c.pd.high")[2,4],2),
                       round(linearHypothesis(c2.e.reg, "c.sfl.high = c.pd.high")[2,4],2)))

myxtable <- xtable(table)

align(myxtable) <- c("l","p{5cm}",rep("p{1.7cm}",6))

caption(myxtable) <- "Potential Spillovers to Women in the Control Group"

label(myxtable) <- "spillovers"

print(myxtable, floating = TRUE, caption.placement="top",sanitize.text.function = identity,
      include.colnames=F,include.rownames=F,table.placement="H",
      hline.after=NULL,
      add.to.row=list(
        pos = list(0,0,1,11,nrow(table)),
        command = c(
          paste0("\\toprule"),
          " & \\multicolumn{2}{c}{Primary Psychosocial Index} & \\multicolumn{2}{c}{Secondary Psychosocial Index} & \\multicolumn{2}{c}{Economic Index}\\\\",
          " \\midrule ",
          " \\midrule ",
          "\\bottomrule \\multicolumn{7}{l}{\\multirow{2}{22.5cm}{\\small Notes: Spillover effects to women in the control group on index variables associated with each of our three families of outcomes. Columns (1), (3), and (5) report effects for each treated woman in the cooperative of a control woman while columns (2), (4), and (6) report effects for control women who have an above-median number of treated women in their cooperative (note that the median is 3 for both the Storytelling and Soft Skills workshops). We report robust standard errors in parentheses. The bottom row shows the p-values associated with hypothesis tests of equality between the estimated average effect for Storytelling (S) versus Soft Skills (SS).}} \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ "
        )),
      type = "latex",file="spillovers.tex")


##########################################################################################
#Online Appendix Figures and Tables  
##########################################################################################
####Figure OA3 - Correlation between MUE and income - as expected, significant and negative#### 
mue.inc.bl <- felm(bl.mue ~ bl.income|0|0, data = full)
summary(mue.inc.bl, robust = TRUE)

mue.inc.fu <- felm(fu.mue ~ fu.income + SFL + PD|0|0, data = full)
summary(mue.inc.fu, robust = TRUE)

mue.wrk.bl <- felm(bl.mue ~ bl.work_lastmonth|0|0, data = full)
summary(mue.wrk.bl, robust = TRUE)

mue.wrk.fu <- felm(fu.mue ~ fu.work_lastmonth + SFL + PD|0|0, data = full)
summary(mue.wrk.fu, robust = TRUE)

#Make a scatterplot of baseline values of log income and MUE
mue.inc.plot <- ggplot(data = full, aes(x = ihs.bl, y = bl.mue)) + geom_point() + 
  geom_smooth(method = "lm", color = "black") + 
  labs(x = "IHS(Cash Earnings)", 
       y = "MUE") + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), 
                     axis.text.x = element_text(angle = 45, hjust = 1)) 

ggsave(filename="Figures/MUE_Inc_Corr.eps", plot=mue.inc.plot, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)

#####Figure OA4 - Treatment effect plots for primary psycho-social outcomes, robustness to no BL cov#####
sa.pd.r <- as.data.frame(cbind(sa.reg.ew2$coefficients, sa.reg.ew2$rse))
sa.pd.r <- cbind(sa.pd.r, rownames(sa.pd.r))
colnames(sa.pd.r) <- c("ATE","RSE","Tmt")
sa.pd.r <- mutate(sa.pd.r, dep.var = "Self-Advocacy Score")

l.pd.r <- as.data.frame(cbind(nlead.reg2$coefficients, nlead.reg2$rse))
l.pd.r <- cbind(l.pd.r, rownames(l.pd.r))
colnames(l.pd.r) <- c("ATE","RSE","Tmt")
l.pd.r <- mutate(l.pd.r, dep.var = "No. Leadership Positions")

lind.pd.r <- as.data.frame(cbind(lead.ind.reg2$coefficients, lead.ind.reg2$rse))
lind.pd.r <- cbind(lind.pd.r, rownames(lind.pd.r))
colnames(lind.pd.r) <- c("ATE","RSE","Tmt")
lind.pd.r <- mutate(lind.pd.r, dep.var = "Any Leadership Position")

cc.pd.r <- as.data.frame(cbind(con.reg.ew2$coefficients, con.reg.ew2$rse))
cc.pd.r <- cbind(cc.pd.r, rownames(cc.pd.r))
colnames(cc.pd.r) <- c("ATE","RSE","Tmt")
cc.pd.r <- mutate(cc.pd.r, dep.var = "Connectedness to Co-op")

chv.pd.r <- as.data.frame(cbind(conhv.reg.ew2$coefficients, conhv.reg.ew2$rse))
chv.pd.r <- cbind(chv.pd.r, rownames(chv.pd.r))
colnames(chv.pd.r) <- c("ATE","RSE","Tmt")
chv.pd.r <- mutate(chv.pd.r, dep.var = "Connectedness to Home/Village")

ps.pd.r <- as.data.frame(cbind(ps.reg.ew2$coefficients, ps.reg.ew2$rse))
ps.pd.r <- cbind(ps.pd.r, rownames(ps.pd.r))
colnames(ps.pd.r) <- c("ATE","RSE","Tmt")
ps.pd.r <- mutate(ps.pd.r, dep.var = "Problem Solving Score")

com.pd.r <- as.data.frame(cbind(com.reg.ew2$coefficients, com.reg.ew2$rse))
com.pd.r <- cbind(com.pd.r, rownames(com.pd.r))
colnames(com.pd.r) <- c("ATE","RSE","Tmt")
com.pd.r <- mutate(com.pd.r, dep.var = "Communication Skills Score")

pps.pd.r <- as.data.frame(cbind(pps.reg2$coefficients, pps.reg2$rse))
pps.pd.r <- cbind(pps.pd.r, rownames(pps.pd.r))
colnames(pps.pd.r) <- c("ATE","RSE","Tmt")
pps.pd.r <- mutate(pps.pd.r, dep.var = "PPS Index")

pps.r.plot.data <- rbind(sa.pd.r, l.pd.r, lind.pd.r, cc.pd.r, chv.pd.r, ps.pd.r, com.pd.r, pps.pd.r) 
pps.r.plot.data <- filter(pps.r.plot.data, Tmt=="SFL"|Tmt=="PD") %>%
  mutate(SD = if_else(dep.var=="Self-Advocacy Score", sa.sd, 
                      if_else(dep.var=="Connectedness to Home/Village", chv.sd,
                              if_else(dep.var=="Problem Solving Score", ps.sd, 
                                      if_else(dep.var=="Communication Skills Score", com.sd,
                                              if_else(dep.var=="Connectedness to Co-op",cc.sd,
                                                      if_else(dep.var=="No. Leadership Positions",l.sd,
                                                              if_else(dep.var=="Any Leadership Position", lind.sd,
                                                                      if_else(dep.var=="PPS Index",pps.sd, as.double(NA))))))))),
         ATE.std = ATE/SD, RSE.std = RSE/SD, CI.high = ATE.std + 1.96*RSE.std,
         CI.low = ATE.std - 1.96*RSE.std,
         CI90.high = ATE.std + 1.65*RSE.std, CI90.low = ATE.std - 1.65*RSE.std)

pps.r.plot.data$dep.var[which(pps.r.plot.data$dep.var=="PPS Index")] <- " PPS Index"
pps.r.plot.data <- mutate(pps.r.plot.data, Treatment = if_else(Tmt=="PD","Soft Skills","Storytelling"))

pps.r.plot <- ggplot(data = pps.r.plot.data, aes(y = dep.var, x = ATE.std, group = Treatment,
                                                 fill = Treatment, color = Treatment, 
                                                 xmin = CI.low, xmax = CI.high)) + 
  geom_point(position = position_dodge(0.5)) + 
  geom_errorbar(width = 0.1, position = position_dodge(0.5)) + 
  geom_vline(xintercept = 0, color = "black") + 
  labs(x = "Treatment Effect (Standard Deviations)", 
       y = "Outcome") + 
  scale_fill_brewer(palette="Paired", direction = -1, aesthetics = c("colour","fill")) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), 
                     axis.text.x = element_text(angle = 45, hjust = 1)) 

ggsave(filename="Figures/PPS_robust.eps", plot=pps.r.plot, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)

#####Figure OA5 - Treatment effect plots for secondary psycho-social outcomes - robustness to no BL cov#####
pss.pd.r <- as.data.frame(cbind(ssp.reg.2$coefficients, ssp.reg.2$rse))
pss.pd.r <- cbind(pss.pd.r, rownames(pss.pd.r))
colnames(pss.pd.r) <- c("ATE","RSE","Tmt")
pss.pd.r <- mutate(pss.pd.r, dep.var = "SSS - Past")

ss.pd.r <- as.data.frame(cbind(ss.reg.2$coefficients, ss.reg.2$rse))
ss.pd.r <- cbind(ss.pd.r, rownames(ss.pd.r))
colnames(ss.pd.r) <- c("ATE","RSE","Tmt")
ss.pd.r <- mutate(ss.pd.r, dep.var = "SSS - Current")

fss.pd.r <- as.data.frame(cbind(ssf.reg.2$coefficients, ssf.reg.2$rse))
fss.pd.r <- cbind(fss.pd.r, rownames(fss.pd.r))
colnames(fss.pd.r) <- c("ATE","RSE","Tmt")
fss.pd.r <- mutate(fss.pd.r, dep.var = "SSS - Future")

conw.pd.r <- as.data.frame(cbind(tg.reg2$coefficients, tg.reg2$rse))
conw.pd.r <- cbind(conw.pd.r, rownames(conw.pd.r))
colnames(conw.pd.r) <- c("ATE","RSE","Tmt")
conw.pd.r <- mutate(conw.pd.r, dep.var = "Connectedness to Workshop Group")

v.pd.r <- as.data.frame(cbind(val.reg.ew2$coefficients, val.reg.ew2$rse))
v.pd.r <- cbind(v.pd.r, rownames(v.pd.r))
colnames(v.pd.r) <- c("ATE","RSE","Tmt")
v.pd.r <- mutate(v.pd.r, dep.var = "Self-Value Score")

r.pd.r <- as.data.frame(cbind(res.reg.ew2$coefficients, res.reg.ew2$rse))
r.pd.r <- cbind(r.pd.r, rownames(r.pd.r))
colnames(r.pd.r) <- c("ATE","RSE","Tmt")
r.pd.r <- mutate(r.pd.r, dep.var = "Resilience Score")

idl.pd.r <- as.data.frame(cbind(lead.reg2$coefficients, lead.reg2$rse))
idl.pd.r <- cbind(idl.pd.r, rownames(idl.pd.r))
colnames(idl.pd.r) <- c("ATE","RSE","Tmt")
idl.pd.r <- mutate(idl.pd.r, dep.var = "ID as Leader Score")

rm.pd.r <- as.data.frame(cbind(rm.reg.2$coefficients, rm.reg.2$rse))
rm.pd.r <- cbind(rm.pd.r, rownames(rm.pd.r))
colnames(rm.pd.r) <- c("ATE","RSE","Tmt")
rm.pd.r <- mutate(rm.pd.r, dep.var = "Has Role Model")

rma.pd.r <- as.data.frame(cbind(rmasp.reg.2$coefficients, rmasp.reg.2$rse))
rma.pd.r <- cbind(rma.pd.r, rownames(rma.pd.r))
colnames(rma.pd.r) <- c("ATE","RSE","Tmt")
rma.pd.r <- mutate(rma.pd.r, dep.var = "Aspires to be Role Model")

rmt.pd.r <- as.data.frame(cbind(rmtime.reg.2$coefficients, rmtime.reg.2$rse))
rmt.pd.r <- cbind(rmt.pd.r, rownames(rmt.pd.r))
colnames(rmt.pd.r) <- c("ATE","RSE","Tmt")
rmt.pd.r <- mutate(rmt.pd.r, dep.var = "Time to be Role Model")

g.pd.r <- as.data.frame(cbind(goal.reg.2$coefficients, goal.reg.2$rse))
g.pd.r <- cbind(g.pd.r, rownames(g.pd.r))
colnames(g.pd.r) <- c("ATE","RSE","Tmt")
g.pd.r <- mutate(g.pd.r, dep.var = "Has a Goal")

gc.pd.r <- as.data.frame(cbind(com.goal.reg$coefficients, com.goal.reg$rse))
gc.pd.r <- cbind(gc.pd.r, rownames(gc.pd.r))
colnames(gc.pd.r) <- c("ATE","RSE","Tmt")
gc.pd.r <- mutate(gc.pd.r, dep.var = "Achieved Last Goal")

sps.pd.r <- as.data.frame(cbind(sps.reg2$coefficients, sps.reg2$rse))
sps.pd.r <- cbind(sps.pd.r, rownames(sps.pd.r))
colnames(sps.pd.r) <- c("ATE","RSE","Tmt")
sps.pd.r <- mutate(sps.pd.r, dep.var = "SPS Index")

sps.r.plot.data <- rbind(pss.pd.r, ss.pd.r, fss.pd.r, conw.pd.r,v.pd.r, r.pd.r, idl.pd.r, rm.pd.r, 
                         rma.pd.r, rmt.pd.r, g.pd.r, gc.pd.r, sps.pd.r) 
sps.r.plot.data <- filter(sps.r.plot.data, Tmt=="SFL"|Tmt=="PD") %>%
  mutate(SD = if_else(dep.var=="SSS - Past", pss.sd, 
                      if_else(dep.var=="SSS - Current", ss.sd,
                              if_else(dep.var=="SSS - Future", fss.sd, 
                                      if_else(dep.var=="Connectedness to Workshop Group", conw.sd,
                                              if_else(dep.var=="Self-Value Score",v.sd,
                                                      if_else(dep.var=="Resilience Score",r.sd,
                                                              if_else(dep.var=="ID as Leader Score", idl.sd,
                                                                      if_else(dep.var=="Has Role Model", rm.sd,
                                                                              if_else(dep.var=="Aspires to be Role Model", rma.sd,
                                                                                      if_else(dep.var=="Time to be Role Model", rmt.sd,
                                                                                              if_else(dep.var=="Has a Goal", g.sd,
                                                                                                      if_else(dep.var=="Achieved Last Goal", gc.sd,
                                                                                                              if_else(dep.var=="SPS Index", sps.sd,
                                                                                                                      as.double(NA)))))))))))))),
         ATE.std = ATE/SD, RSE.std = RSE/SD, CI.high = ATE.std + 1.96*RSE.std,
         CI.low = ATE.std - 1.96*RSE.std,
         CI90.high = ATE.std + 1.65*RSE.std, CI90.low = ATE.std - 1.65*RSE.std)

sps.r.plot.data$dep.var <- factor(sps.r.plot.data$dep.var, 
                                  levels = c("SPS Index","Achieved Last Goal","Has a Goal",
                                             "Time to be Role Model","Aspires to be Role Model",
                                             "Has Role Model","ID as Leader Score","Resilience Score",
                                             "Self-Value Score", "Connectedness to Workshop Group", 
                                             "SSS - Future","SSS - Past","SSS - Current"))
sps.r.plot.data <- mutate(sps.r.plot.data, Treatment = if_else(Tmt=="PD","Soft Skills","Storytelling"))

sps.r.plot.1 <- ggplot(data = sps.r.plot.data, 
                       aes(y = dep.var, x = ATE.std, group = Treatment,
                           fill = Treatment, color = Treatment, 
                           xmin = CI.low, xmax = CI.high)) + 
  geom_point(position = position_dodge(0.5)) + 
  geom_errorbar(width = 0.1, position = position_dodge(0.5)) + 
  geom_vline(xintercept = 0, color = "black") + 
  labs(x = "Treatment Effect (Standard Deviations)", 
       y = "Outcome") + 
  scale_fill_brewer(palette="Paired", direction = -1, aesthetics = c("colour","fill")) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), 
                     axis.text.x = element_text(angle = 45, hjust = 1)) 

ggsave(filename="Figures/SPS_robust.eps", plot=sps.r.plot.1, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)

#####Figure OA6 - Treatment effect plots for economic outcomes - robustness to no BL cov#####
inc.pd.r <- as.data.frame(cbind(inc2$coefficients, inc2$rse))
inc.pd.r <- cbind(inc.pd.r, rownames(inc.pd.r))
colnames(inc.pd.r) <- c("ATE","RSE","Tmt")
inc.pd.r <- mutate(inc.pd.r, dep.var = "Cash Earnings")

inci.pd.r <- as.data.frame(cbind(work.reg2$coefficients, work.reg2$rse))
inci.pd.r <- cbind(inci.pd.r, rownames(inci.pd.r))
colnames(inci.pd.r) <- c("ATE","RSE","Tmt")
inci.pd.r <- mutate(inci.pd.r, dep.var = "Earn Any Cash")

mue.pd.r <- as.data.frame(cbind(mue.reg2$coefficients, mue.reg2$rse))
mue.pd.r <- cbind(mue.pd.r, rownames(mue.pd.r))
colnames(mue.pd.r) <- c("ATE","RSE","Tmt")
mue.pd.r <- mutate(mue.pd.r, dep.var = "Marginal Utility of Expenditure")

exp.pd.r <- as.data.frame(cbind(exp.reg2$coefficients, exp.reg2$rse))
exp.pd.r <- cbind(exp.pd.r, rownames(exp.pd.r))
colnames(exp.pd.r) <- c("ATE","RSE","Tmt")
exp.pd.r <- mutate(exp.pd.r, dep.var = "Log(Consumption Expenditures per Capita)")

e.pd.r <- as.data.frame(cbind(e.reg2$coefficients, e.reg2$rse))
e.pd.r <- cbind(e.pd.r, rownames(e.pd.r))
colnames(e.pd.r) <- c("ATE","RSE","Tmt")
e.pd.r <- mutate(e.pd.r, dep.var = "Econ Index")

e.r.plot.data <- rbind(inc.pd.r, inci.pd.r, mue.pd.r, exp.pd.r, e.pd.r) 
e.r.plot.data <- filter(e.r.plot.data, Tmt=="SFL"|Tmt=="PD") %>%
  mutate(SD = if_else(dep.var=="Cash Earnings", inc.sd, 
                      if_else(dep.var=="Earn Any Cash", inci.sd,
                              if_else(dep.var=="Marginal Utility of Expenditure", mue.sd,
                                      if_else(dep.var=="Log(Consumption Expenditures per Capita)",exp.sd, 
                                              if_else(dep.var=="Econ Index", e.sd, as.double(NA)))))),
         ATE.std = ATE/SD, RSE.std = RSE/SD, CI.high = ATE.std + 1.96*RSE.std,
         CI.low = ATE.std - 1.96*RSE.std)

e.r.plot.data$dep.var[which(e.r.plot.data$dep.var=="Econ Index")] <- " Econ Index"
e.r.plot.data <- mutate(e.r.plot.data, Treatment = if_else(Tmt=="PD", "Soft Skills","Storytelling"))

e.r.plot <- ggplot(data = e.r.plot.data, aes(y = dep.var, x = ATE.std, group = Treatment,
                                             fill = Treatment, color = Treatment, 
                                             xmin = CI.low, xmax = CI.high)) + 
  geom_point(position = position_dodge(0.5)) + 
  geom_errorbar(width = 0.1, position = position_dodge(0.5)) + 
  geom_vline(xintercept = 0, color = "black") + 
  labs(x = "Treatment Effect (Standard Deviations)", 
       y = "Outcome") + 
  scale_fill_brewer(palette="Paired", direction = -1, aesthetics = c("colour","fill")) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), 
                     axis.text.x = element_text(angle = 45, hjust = 1)) 

ggsave(filename="Figures/Econ_robust.eps", plot=e.r.plot, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)




#####Table OA1 - Balance on secondary psychosocial outcomes#####
table <- matrix(NA, nrow=24, ncol=5)

table[1,] <- c("","Cash Control", "Soft Skills", "Storytelling","p-value")
table[2,] <- unlist(c("SSS - Current", balance.means[10,1:3],round(unlist(summary(bal.ss, robust = TRUE)[11])[4],2)))
table[3,] <- c("",paste0("(",balance.sds[10,1],")"),paste0("(",balance.sds[10,2],")"),paste0("(",balance.sds[10,3],")"),"")
table[4,] <- unlist(c("SSS - Past", balance.means[11,1:3],round(unlist(summary(bal.pss, robust = TRUE)[11])[4],2)))
table[5,] <- c("",paste0("(",balance.sds[11,1],")"),paste0("(",balance.sds[11,2],")"),paste0("(",balance.sds[11,3],")"),"")
table[6,] <- unlist(c("SSS - Future", balance.means[12,1:3],round(unlist(summary(bal.fss, robust = TRUE)[11])[4],2)))
table[7,] <- c("",paste0("(",balance.sds[12,1],")"),paste0("(",balance.sds[12,2],")"),paste0("(",balance.sds[12,3],")"),"")
table[8,] <- unlist(c("Self-Value", balance.means[13,1:3],round(unlist(summary(bal.v, robust = TRUE)[11])[4],2)))
table[9,] <- c("",paste0("(",balance.sds[13,1],")"),paste0("(",balance.sds[13,2],")"),paste0("(",balance.sds[13,3],")"),"")
table[10,] <- unlist(c("Resilience", balance.means[14,1:3],round(unlist(summary(bal.r, robust = TRUE)[11])[4],2)))
table[11,] <- c("",paste0("(",balance.sds[14,1],")"),paste0("(",balance.sds[14,2],")"),paste0("(",balance.sds[14,3],")"),"")
table[12,] <- unlist(c("Identify as Leader", balance.means[15,1:3],round(unlist(summary(bal.lead, robust = TRUE)[11])[4],2)))
table[13,] <- c("",paste0("(",balance.sds[15,1],")"),paste0("(",balance.sds[15,2],")"),paste0("(",balance.sds[15,3],")"),"")
table[14,] <- unlist(c("Connectedness to Workshop Group", balance.means[16,1:3],round(unlist(summary(bal.tg, robust = TRUE)[11])[4],2)))
table[15,] <- c("",paste0("(",balance.sds[16,1],")"),paste0("(",balance.sds[16,2],")"),paste0("(",balance.sds[16,3],")"),"")
table[16,] <- unlist(c("Has Role Model", balance.means[17,1:3],round(unlist(summary(bal.rm, robust = TRUE)[11])[4],2)))
table[17,] <- c("",paste0("(",balance.sds[17,1],")"),paste0("(",balance.sds[17,2],")"),paste0("(",balance.sds[17,3],")"),"")
table[18,] <- unlist(c("Aspires to be Role Model", balance.means[18,1:3],round(unlist(summary(bal.rmasp, robust = TRUE)[11])[4],2)))
table[19,] <- c("",paste0("(",balance.sds[18,1],")"),paste0("(",balance.sds[18,2],")"),paste0("(",balance.sds[18,3],")"),"")
table[20,] <- unlist(c("Time to be Role Model", balance.means[19,1:3],round(unlist(summary(bal.rmtime, robust = TRUE)[11])[4],2)))
table[21,] <- c("",paste0("(",balance.sds[19,1],")"),paste0("(",balance.sds[19,2],")"),paste0("(",balance.sds[19,3],")"),"")
table[22,] <- unlist(c("Has a goal", balance.means[20,1:3],round(unlist(summary(bal.goal, robust = TRUE)[11])[4],2)))
table[23,] <- c("",paste0("(",balance.sds[20,1],")"),paste0("(",balance.sds[20,2],")"),paste0("(",balance.sds[20,3],")"),"")
table[24,] <- unlist(c("N",length(bl$resp_id[which(bl$workshop_name=="Control")]),length(bl$resp_id[which(bl$workshop_name=="Professional Development")]),
                       length(bl$resp_id[which(bl$workshop_name=="Storytelling for Leadership")]),""))


myxtable <- xtable(table)

align(myxtable) <- c("l",rep("c",5))

caption(myxtable) <- "Balance: Pre-Registered Secondary Psychosocial Outcomes"

label(myxtable) <- "balance2"

print(myxtable, floating = TRUE, caption.placement="top",sanitize.text.function = identity,
      include.colnames=F,include.rownames=F,table.placement="H",
      hline.after=NULL,
      add.to.row=list(
        pos = list(0,1,nrow(table)-1,nrow(table)),
        command = c(
          paste0("\\toprule"),
          " \\midrule ",
          " \\midrule ",
          "\\bottomrule \\multicolumn{5}{l}{\\multirow{2}{14cm}{\\small Notes: Mean baseline covariates by treatment group. Standard deviations are in parentheses. Column 4 reports p-values associated with F-tests of joint equality between the three groups.}} \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ "
        )),
      type = "latex",file="balance2.tex")

#####Figure OA7 - Post workshop treatment effects#####
#Connection to training group
full <- mutate(full, pw.tg = if_else(workshop_name=="Control",bl.tg, pw.tg))

tg.reg2.test <- felm(pw.tg ~ SFL + PD|0|0, data = full)
summary(tg.reg2.test, robust = TRUE)

#Self-value
full <- mutate(full, pw.vew = if_else(workshop_name=="Control", bl.vew, pw.vew))

val.reg.ew2.test <- felm(pw.vew ~ SFL + PD |0|0, data = full)
summary(val.reg.ew2.test, robust = TRUE)

#Resilience
full <- mutate(full, pw.rew = if_else(workshop_name=="Control",bl.rew, pw.rew))

res.reg.ew2.test <- felm(pw.rew ~ SFL + PD|0|0, data = full)
summary(res.reg.ew2.test, robust = TRUE)

#Problem solving
full <- mutate(full, pw.psew = if_else(workshop_name=="Control",bl.psew,pw.psew))

ps.reg.ew2.test <- felm(pw.psew ~ SFL + PD|0|0, data = full)
summary(ps.reg.ew2.test, robust = TRUE)

#ID as leader
full <- mutate(full, pw.lead = if_else(workshop_name=="Control",bl.lead, pw.lead))

lead.reg.test <- felm(pw.lead ~ SFL + PD|0|0, data = full)
summary(lead.reg.test, robust = TRUE)

#Have a role model
full <- mutate(full, pw.rm = if_else(workshop_name=="Control", bl.rm, pw.rm))

rm.reg.1.2.test <- felm(pw.rm ~ SFL + PD |0|0, data = full)
summary(rm.reg.1.2.test, robust = TRUE)

#Aspire to be as successful as role model
full <- mutate(full, pw.rmasp = if_else(workshop_name=="Control", bl.rmasp, pw.rmasp))

rm.reg.2.2.test <- felm(pw.rmasp ~ SFL + PD|0|0, data = full)
summary(rm.reg.2.2.test, robust = TRUE)

#Time to achieve role model
full <- mutate(full, pw.rmtime = if_else(workshop_name=="Control", bl.rmtime, pw.rmtime))

rm.reg.3.2.test <- felm(pw.rmtime ~ SFL + PD |0|0, data = full)
summary(rm.reg.3.2.test, robust = TRUE)

tg.sd.pw <- sd(full$bl.tg[which(full$workshop_name=="Control")], na.rm=T)
tg.pd.pw <- as.data.frame(cbind(tg.reg2.test$coefficients, tg.reg2.test$rse))
tg.pd.pw <- cbind(tg.pd.pw, rownames(tg.pd.pw))
colnames(tg.pd.pw) <- c("ATE","RSE","Tmt")
tg.pd.pw <- mutate(tg.pd.pw, dep.var = "Connectedness to Training Group")

v.sd.pw <- sd(full$bl.vew[which(full$workshop_name=="Control")], na.rm=T)
v.pd.pw <- as.data.frame(cbind(val.reg.ew2.test$coefficients, val.reg.ew2.test$rse))
v.pd.pw <- cbind(v.pd.pw, rownames(v.pd.pw))
colnames(v.pd.pw) <- c("ATE","RSE","Tmt")
v.pd.pw <- mutate(v.pd.pw, dep.var = "Self-value")

r.sd.pw <- sd(full$bl.rew[which(full$workshop_name=="Control")], na.rm=T)
r.pd.pw <- as.data.frame(cbind(res.reg.ew2.test$coefficients, res.reg.ew2.test$rse))
r.pd.pw <- cbind(r.pd.pw, rownames(r.pd.pw))
colnames(r.pd.pw) <- c("ATE","RSE","Tmt")
r.pd.pw <- mutate(r.pd.pw, dep.var = "Resilience")

ps.sd.pw <- sd(full$bl.psew[which(full$workshop_name=="Control")], na.rm=T)
ps.pd.pw <- as.data.frame(cbind(ps.reg.ew2.test$coefficients, ps.reg.ew2.test$rse))
ps.pd.pw <- cbind(ps.pd.pw, rownames(ps.pd.pw))
colnames(ps.pd.pw) <- c("ATE","RSE","Tmt")
ps.pd.pw <- mutate(ps.pd.pw, dep.var = "Problem Solving")

l.sd.pw <- sd(full$bl.lead[which(full$workshop_name=="Control")], na.rm=T)
l.pd.pw <- as.data.frame(cbind(lead.reg.test$coefficients, lead.reg.test$rse))
l.pd.pw <- cbind(l.pd.pw, rownames(l.pd.pw))
colnames(l.pd.pw) <- c("ATE","RSE","Tmt")
l.pd.pw <- mutate(l.pd.pw, dep.var = "Identify as Leader")

rm.sd.pw <- sqrt(mean(full$bl.rm[which(full$workshop_name=="Control")], na.rm=T)*(1-mean(full$bl.rm[which(full$workshop_name=="Control")], na.rm=T)))
rm.pd.pw <- as.data.frame(cbind(rm.reg.1.2.test$coefficients, rm.reg.1.2.test$rse))
rm.pd.pw <- cbind(rm.pd.pw, rownames(rm.pd.pw))
colnames(rm.pd.pw) <- c("ATE","RSE","Tmt")
rm.pd.pw <- mutate(rm.pd.pw, dep.var = "Has a Role Model")

rm2.sd.pw <- sqrt(mean(full$bl.rmasp[which(full$workshop_name=="Control")], na.rm=T)*(1-mean(full$bl.rmasp[which(full$workshop_name=="Control")], na.rm=T)))
rm2.pd.pw <- as.data.frame(cbind(rm.reg.2.2.test$coefficients, rm.reg.2.2.test$rse))
rm2.pd.pw <- cbind(rm2.pd.pw, rownames(rm2.pd.pw))
colnames(rm2.pd.pw) <- c("ATE","RSE","Tmt")
rm2.pd.pw <- mutate(rm2.pd.pw, dep.var = "Aspires to be Role Model")

rm3.sd.pw <- sd(full$bl.rmtime[which(full$workshop_name=="Control")], na.rm=T)
rm3.pd.pw <- as.data.frame(cbind(rm.reg.3.2.test$coefficients, rm.reg.3.2.test$rse))
rm3.pd.pw <- cbind(rm3.pd.pw, rownames(rm3.pd.pw))
colnames(rm3.pd.pw) <- c("ATE","RSE","Tmt")
rm3.pd.pw <- mutate(rm3.pd.pw, dep.var = "Time to be Role Model")

pw.plot.data <- rbind(tg.pd.pw, v.pd.pw, r.pd.pw, ps.pd.pw, l.pd.pw, rm.pd.pw, rm2.pd.pw, rm3.pd.pw) 
pw.plot.data <- filter(pw.plot.data, Tmt=="SFL"|Tmt=="PD") %>%
  mutate(SD = if_else(dep.var=="Connectedness to Training Group", tg.sd.pw, 
                      if_else(dep.var=="Self-value", v.sd.pw,
                              if_else(dep.var=="Problem Solving", ps.sd.pw, 
                                      if_else(dep.var=="Resilience", r.sd.pw,
                                              if_else(dep.var=="Identify as Leader",l.sd.pw,
                                                      if_else(dep.var=="Has a Role Model",r.sd.pw,
                                                              if_else(dep.var=="Aspires to be Role Model", rm2.sd.pw,
                                                                      if_else(dep.var=="Time to be Role Model", rm3.sd.pw,
                                                                              as.double(NA))))))))),
         ATE.std = ATE/SD, RSE.std = RSE/SD, CI.high = ATE.std + 1.96*RSE.std,
         CI.low = ATE.std - 1.96*RSE.std,
         CI90.high = ATE.std + 1.65*RSE.std, CI90.low = ATE.std - 1.65*RSE.std)
pw.plot.data <- mutate(pw.plot.data, Treatment = if_else(Tmt=="PD","Soft Skills","Storytelling"))

pw.plot <- ggplot(data = pw.plot.data, aes(y = dep.var, x = ATE.std, group = Treatment,
                                           fill = Treatment, color = Treatment, 
                                           xmin = CI.low, xmax = CI.high)) + 
  geom_point(position = position_dodge(0.5)) + 
  geom_errorbar(width = 0.1, position = position_dodge(0.5)) + 
  geom_vline(xintercept = 0, color = "black") + 
  labs(x = "Treatment Effect (Standard Deviations)", 
       y = "Outcome") + 
  scale_fill_brewer(palette="Paired", direction = -1, aesthetics = c("colour","fill")) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), 
                     axis.text.x = element_text(angle = 45, hjust = 1)) 

ggsave(filename="Figures/PW.eps", plot=pw.plot, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)



####Table OA2 - Heterogeneity by above versus below baseline PPS index ####
table <- matrix(NA, nrow=17, ncol=4)

table[1,] <- c("","Primary Psychosocial Index","Secondary Psychosocial Index"," Economic Index")
table[2,] <- unlist(c("Storytelling", round(het.pps$coefficients[2],2), round(het.sps$coefficients[2],2),
                      round(het.e$coefficients[2],2)))
table[3,] <- c("",paste0("(",round(het.pps$rse[2],2),")"), paste0("(",round(het.sps$rse[2],2),")"),
               paste0("(",round(het.e$rse[2],2),")"))
table[4,] <- unlist(c("Soft Skills", round(het.pps$coefficients[3],2), round(het.sps$coefficients[3],2),
                      round(het.e$coefficients[3],2)))
table[5,] <- c("",paste0("(",round(het.pps$rse[3],2),")"), paste0("(",round(het.sps$rse[3],2),")"),
               paste0("(",round(het.e$rse[3],2),")"))
table[6,] <- unlist(c("Baseline Level", round(het.pps$coefficients[4],2), round(het.sps$coefficients[4],2),
                      round(het.e$coefficients[4],2)))
table[7,] <- c("",paste0("(",round(het.pps$rse[4],2),")"), paste0("(",round(het.sps$rse[4],2),")"),
               paste0("(",round(het.e$rse[4],2),")"))
table[8,] <- unlist(c("Above-Median PPS", round(het.pps$coefficients[5],2), round(het.sps$coefficients[5],2),
                      round(het.e$coefficients[5],2)))
table[9,] <- c("",paste0("(",round(het.pps$rse[5],2),")"), paste0("(",round(het.sps$rse[5],2),")"),
               paste0("(",round(het.e$rse[5],2),")"))
table[10,] <- unlist(c("Storytelling x", round(het.pps$coefficients[6],2), round(het.sps$coefficients[6],2),
                       round(het.e$coefficients[6],2)))
table[11,] <- c("Above-Median PPS",paste0("(",round(het.pps$rse[6],2),")"), paste0("(",round(het.sps$rse[6],2),")"),
                paste0("(",round(het.e$rse[6],2),")"))
table[12,] <- unlist(c("Soft Skills x", round(het.pps$coefficients[7],2), round(het.sps$coefficients[7],2),
                       round(het.e$coefficients[7],2)))
table[13,] <- c("Above-Median PPS",paste0("(",round(het.pps$rse[7],2),")"), paste0("(",round(het.sps$rse[7],2),")"),
                paste0("(",round(het.e$rse[7],2),")"))
table[14,] <- unlist(c("N",het.pps$N, het.sps$N, het.e$N))
table[15,] <- c("Adj. R2",round(as.numeric(summary(het.pps)[15]),2), round(as.numeric(summary(het.sps)[15]),2),
                round(as.numeric(summary(het.e)[15]),2))
table[16,] <- unlist(c("S=SS (p-value)", round(linearHypothesis(het.pps, "SFL = PD")[2,4],2), round(linearHypothesis(het.sps, "SFL = PD")[2,4],2),
                       round(linearHypothesis(het.e, "SFL = PD")[2,4],2)))
table[17,] <- unlist(c("S x A-M PPS=SS x A-M PPS (p-value)", round(linearHypothesis(het.pps, "SFL:high.pps = PD:high.pps")[2,4],2), 
                       round(linearHypothesis(het.sps, "SFL:high.pps = PD:high.pps")[2,4],2),
                       round(linearHypothesis(het.e, "SFL:high.pps = PD:high.pps")[2,4],2)))

myxtable <- xtable(table)

align(myxtable) <- c("l","p{5cm}",rep("p{3cm}",3))

caption(myxtable) <- "Heterogeneity by Above-Median Baseline Primary Psychosocial Index Values"

label(myxtable) <- "ppshet"

print(myxtable, floating = TRUE, caption.placement="top",sanitize.text.function = identity,
      include.colnames=F,include.rownames=F,table.placement="H",
      hline.after=NULL,
      add.to.row=list(
        pos = list(0,1,13,nrow(table)),
        command = c(
          paste0("\\toprule"),
          " \\midrule ",
          " \\midrule ",
          "\\bottomrule \\multicolumn{4}{l}{\\multirow{2}{16.5cm}{\\small Notes: Average treatment effects on index variables associated with each of our three families of outcomes, heterogeneously by below- versus above-median values on the index of primary psychosocial outcomes at baseline. We report robust standard errors in parentheses. The bottom two rows show the p-values associated with hypothesis tests of equality between the estimated average treatment effect for Storytelling (S) and Soft Skills (SS), heterogeneously.}} \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ "
        )),
      type = "latex",file="ppshet.tex")

#####Table OA3 - Heterogeneity by above versus below baseline SPS index#####
table <- matrix(NA, nrow=17, ncol=4)

table[1,] <- c("","Primary Psychosocial Index","Secondary Psychosocial Index"," Economic Index")
table[2,] <- unlist(c("Storytelling", round(het2.pps$coefficients[2],2), round(het2.sps$coefficients[2],2),
                      round(het2.e$coefficients[2],2)))
table[3,] <- c("",paste0("(",round(het2.pps$rse[2],2),")"), paste0("(",round(het2.sps$rse[2],2),")"),
               paste0("(",round(het2.e$rse[2],2),")"))
table[4,] <- unlist(c("Soft Skills", round(het2.pps$coefficients[3],2), round(het2.sps$coefficients[3],2),
                      round(het2.e$coefficients[3],2)))
table[5,] <- c("",paste0("(",round(het2.pps$rse[3],2),")"), paste0("(",round(het2.sps$rse[3],2),")"),
               paste0("(",round(het2.e$rse[3],2),")"))
table[6,] <- unlist(c("Baseline Level", round(het2.pps$coefficients[4],2), round(het2.sps$coefficients[4],2),
                      round(het2.e$coefficients[4],2)))
table[7,] <- c("",paste0("(",round(het2.pps$rse[4],2),")"), paste0("(",round(het2.sps$rse[4],2),")"),
               paste0("(",round(het2.e$rse[4],2),")"))
table[8,] <- unlist(c("Above-Median SPS", round(het2.pps$coefficients[5],2), round(het2.sps$coefficients[5],2),
                      round(het2.e$coefficients[5],2)))
table[9,] <- c("",paste0("(",round(het2.pps$rse[5],2),")"), paste0("(",round(het2.sps$rse[5],2),")"),
               paste0("(",round(het2.e$rse[5],2),")"))
table[10,] <- unlist(c("Storytelling x", round(het2.pps$coefficients[6],2), round(het2.sps$coefficients[6],2),
                       round(het2.e$coefficients[6],2)))
table[11,] <- c("Above-Median SPS",paste0("(",round(het2.pps$rse[6],2),")"), paste0("(",round(het2.sps$rse[6],2),")"),
                paste0("(",round(het2.e$rse[6],2),")"))
table[12,] <- unlist(c("Soft Skills x", round(het2.pps$coefficients[7],2), round(het2.sps$coefficients[7],2),
                       round(het2.e$coefficients[7],2)))
table[13,] <- c("Above-Median SPS",paste0("(",round(het2.pps$rse[7],2),")"), paste0("(",round(het2.sps$rse[7],2),")"),
                paste0("(",round(het2.e$rse[7],2),")"))
table[14,] <- unlist(c("N",het2.pps$N, het2.sps$N, het2.e$N))
table[15,] <- c("Adj. R2",round(as.numeric(summary(het2.pps)[15]),2), round(as.numeric(summary(het2.sps)[15]),2),
                round(as.numeric(summary(het2.e)[15]),2))
table[16,] <- unlist(c("S=SS (p-value)", round(linearHypothesis(het2.pps, "SFL = PD")[2,4],2), round(linearHypothesis(het2.sps, "SFL = PD")[2,4],2),
                       round(linearHypothesis(het2.e, "SFL = PD")[2,4],2)))
table[17,] <- unlist(c("S x A-M SPS=SS x A-M SPS (p-value)", round(linearHypothesis(het2.pps, "SFL:high.sps = PD:high.sps")[2,4],2), 
                       round(linearHypothesis(het2.sps, "SFL:high.sps = PD:high.sps")[2,4],2),
                       round(linearHypothesis(het2.e, "SFL:high.sps = PD:high.sps")[2,4],2)))

myxtable <- xtable(table)

align(myxtable) <- c("l","p{5cm}",rep("p{3cm}",3))

caption(myxtable) <- "Heterogeneity by Above-Median Baseline Secondary Psychosocial Index Values"

label(myxtable) <- "spshet"

print(myxtable, floating = TRUE, caption.placement="top",sanitize.text.function = identity,
      include.colnames=F,include.rownames=F,table.placement="H",
      hline.after=NULL,
      add.to.row=list(
        pos = list(0,1,13,nrow(table)),
        command = c(
          paste0("\\toprule"),
          " \\midrule ",
          " \\midrule ",
          "\\bottomrule \\multicolumn{4}{l}{\\multirow{2}{16.5cm}{\\small Notes: Average treatment effects on index variables associated with each of our three families of outcomes, heterogeneously by below- versus above-median values on the index of secondary psychosocial outcomes at baseline. We report robust standard errors in parentheses. The bottom two rows show the p-values associated with hypothesis tests of equality between the estimated average treatment effect for Storytelling (S) versus Soft Skills (SS), heterogeneously.}} \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ "
        )),
      type = "latex",file="spshet.tex")

#####Table OA4 - Heterogeneity by above versus below baseline Econ index######
table <- matrix(NA, nrow=17, ncol=4)

table[1,] <- c("","Primary Psychosocial Index","Secondary Psychosocial Index"," Economic Index")
table[2,] <- unlist(c("Storytelling", round(het3.pps$coefficients[2],2), round(het3.sps$coefficients[2],2),
                      round(het3.e$coefficients[2],2)))
table[3,] <- c("",paste0("(",round(het3.pps$rse[2],2),")"), paste0("(",round(het3.sps$rse[2],2),")"),
               paste0("(",round(het3.e$rse[2],2),")"))
table[4,] <- unlist(c("Soft Skills", round(het3.pps$coefficients[3],2), round(het3.sps$coefficients[3],2),
                      round(het3.e$coefficients[3],2)))
table[5,] <- c("",paste0("(",round(het3.pps$rse[3],2),")"), paste0("(",round(het3.sps$rse[3],2),")"),
               paste0("(",round(het3.e$rse[3],2),")"))
table[6,] <- unlist(c("Baseline Level", round(het3.pps$coefficients[4],2), round(het3.sps$coefficients[4],2),
                      round(het3.e$coefficients[4],2)))
table[7,] <- c("",paste0("(",round(het3.pps$rse[4],2),")"), paste0("(",round(het3.sps$rse[4],2),")"),
               paste0("(",round(het3.e$rse[4],2),")"))
table[8,] <- unlist(c("Above-Median Econ", round(het3.pps$coefficients[5],2), round(het3.sps$coefficients[5],2),
                      round(het3.e$coefficients[5],2)))
table[9,] <- c("",paste0("(",round(het3.pps$rse[5],2),")"), paste0("(",round(het3.sps$rse[5],2),")"),
               paste0("(",round(het3.e$rse[5],2),")"))
table[10,] <- unlist(c("Storytelling x", round(het3.pps$coefficients[6],2), round(het3.sps$coefficients[6],2),
                       round(het3.e$coefficients[6],2)))
table[11,] <- c("Above-Median Econ",paste0("(",round(het3.pps$rse[6],2),")"), paste0("(",round(het3.sps$rse[6],2),")"),
                paste0("(",round(het3.e$rse[6],2),")"))
table[12,] <- unlist(c("Soft Skills x", round(het3.pps$coefficients[7],2), round(het3.sps$coefficients[7],2),
                       round(het3.e$coefficients[7],2)))
table[13,] <- c("Above-Median Econ",paste0("(",round(het3.pps$rse[7],2),")"), paste0("(",round(het3.sps$rse[7],2),")"),
                paste0("(",round(het3.e$rse[7],2),")"))
table[14,] <- unlist(c("N",het3.pps$N, het3.sps$N, het3.e$N))
table[15,] <- c("Adj. R2",round(as.numeric(summary(het3.pps)[15]),2), round(as.numeric(summary(het3.sps)[15]),2),
                round(as.numeric(summary(het3.e)[15]),2))
table[16,] <- unlist(c("S=SS (p-value)", round(linearHypothesis(het3.pps, "SFL = PD")[2,4],2), round(linearHypothesis(het3.sps, "SFL = PD")[2,4],2),
                       round(linearHypothesis(het3.e, "SFL = PD")[2,4],2)))
table[17,] <- unlist(c("S x A-M Econ=SS x A-M Econ (p-value)", round(linearHypothesis(het3.pps, "SFL:high.e = PD:high.e")[2,4],2), 
                       round(linearHypothesis(het3.sps, "SFL:high.e = PD:high.e")[2,4],2),
                       round(linearHypothesis(het3.e, "SFL:high.e = PD:high.e")[2,4],2)))

myxtable <- xtable(table)

align(myxtable) <- c("l","p{5cm}",rep("p{3cm}",3))

caption(myxtable) <- "Heterogeneity by Above-Median Baseline Secondary Economic Index Values"

label(myxtable) <- "ehet"

print(myxtable, floating = TRUE, caption.placement="top",sanitize.text.function = identity,
      include.colnames=F,include.rownames=F,table.placement="H",
      hline.after=NULL,
      add.to.row=list(
        pos = list(0,1,13,nrow(table)),
        command = c(
          paste0("\\toprule"),
          " \\midrule ",
          " \\midrule ",
          "\\bottomrule \\multicolumn{4}{l}{\\multirow{2}{16.5cm}{\\small Notes: Average treatment effects on index variables associated with each of our three families of outcomes, heterogeneously by below- versus above-median values on the index of economic outcomes at baseline. We report robust standard errors in parentheses. The bottom two rows show the p-values associated with hypothesis tests of equality between the estimated average treatment effect for Storytelling (S) versus Soft Skills (SS), heterogeneously.}} \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ "
        )),
      type = "latex",file="ehet.tex")


#####Figure OA8 - Heterogeneity by Age, PPS####
u40.sfl <- as.data.frame(rbind(c(sa.het$coefficients[2], sa.het$rse[2]),
                               c(ps.het$coefficients[2], ps.het$rse[2]),
                               c(nl.het$coefficients[2], nl.het$rse[2]),
                               c(chv.het$coefficients[2], chv.het$rse[2]),
                               c(cc.het$coefficients[2], cc.het$rse[2]),
                               c(com.het$coefficients[2], com.het$rse[2]),
                               c(li.het$coefficients[2], li.het$rse[2]),
                               c(pps.het$coefficients[2], pps.het$rse[2])))

o40.sfl <- as.data.frame(rbind(c(sa.het2$coefficients[2], sa.het2$rse[2]),
                               c(ps.het2$coefficients[2], ps.het2$rse[2]),
                               c(nl.het2$coefficients[2], nl.het2$rse[2]),
                               c(chv.het2$coefficients[2], chv.het2$rse[2]),
                               c(cc.het2$coefficients[2], cc.het2$rse[2]),
                               c(com.het2$coefficients[2], com.het2$rse[2]),
                               c(li.het2$coefficients[2], li.het2$rse[2]),
                               c(pps.het2$coefficients[2], pps.het$rse[2])))

u40.pd <- as.data.frame(rbind(c(sa.het$coefficients[3], sa.het$rse[3]),
                              c(ps.het$coefficients[3], ps.het$rse[3]),
                              c(nl.het$coefficients[3], nl.het$rse[3]),
                              c(chv.het$coefficients[3], chv.het$rse[3]),
                              c(cc.het$coefficients[3], cc.het$rse[3]),
                              c(com.het$coefficients[3], com.het$rse[3]),
                              c(li.het$coefficients[3], li.het$rse[3]),
                              c(pps.het$coefficients[3], pps.het$rse[3])))

o40.pd <- as.data.frame(rbind(c(sa.het2$coefficients[3], sa.het2$rse[3]),
                              c(ps.het2$coefficients[3], ps.het2$rse[3]),
                              c(nl.het2$coefficients[3], nl.het2$rse[3]),
                              c(chv.het2$coefficients[3], chv.het2$rse[3]),
                              c(cc.het2$coefficients[3], cc.het2$rse[3]),
                              c(com.het2$coefficients[3], com.het2$rse[3]),
                              c(li.het2$coefficients[3], li.het2$rse[3]),
                              c(pps.het2$coefficients[3], pps.het2$rse[3])))

cash.u40 <- filter(full, workshop_name=="Control" & under40==1)
cash.o40 <- filter(full, workshop_name=="Control" & over40==1)

u40.sds <- as.data.frame(rbind(c(sd(cash.u40$fu.saew,na.rm=T)),
                               c(sd(cash.u40$fu.psew,na.rm=T)),
                               c(sd(cash.u40$fu.nlead,na.rm=T)),
                               c(sd(cash.u40$fu.conhvew,na.rm=T)),
                               c(sd(cash.u40$fu.conew,na.rm=T)),
                               c(sd(cash.u40$fu.comew,na.rm=T)),
                               c(sqrt(mean(cash.u40$fu.lead.ind,na.rm=T)*(1-mean(cash.u40$fu.lead.ind,na.rm=T)))),
                               c(sd(cash.u40$fu.pps.ind, na.rm=T))))

o40.sds <- as.data.frame(rbind(c(sd(cash.o40$fu.saew,na.rm=T)),
                               c(sd(cash.o40$fu.psew,na.rm=T)),
                               c(sd(cash.o40$fu.nlead,na.rm=T)),
                               c(sd(cash.o40$fu.conhvew,na.rm=T)),
                               c(sd(cash.o40$fu.conew,na.rm=T)),
                               c(sd(cash.o40$fu.comew,na.rm=T)),
                               c(sqrt(mean(cash.o40$fu.lead.ind,na.rm=T)*(1-mean(cash.o40$fu.lead.ind,na.rm=T)))),
                               c(sd(cash.o40$fu.pps.ind, na.rm=T))))

pps.out <- as.data.frame(c("Self-Advocacy Score","Problem Solving Score",
                           "No. Leadership Positions","Connectedness to Home/Village",
                           "Connectedness to Co-op","Communication Skills Score",
                           "Any Leadership Position"," PPS Index"))

u40.sfl <- cbind(u40.sfl, u40.sds, pps.out)
colnames(u40.sfl) <- c("ATE","RSE","Cash.SD","dep.var")
u40.sfl <- mutate(u40.sfl, age.group = "Under 40", tmt="SFL")

o40.sfl <- cbind(o40.sfl, o40.sds, pps.out)
colnames(o40.sfl) <- c("ATE","RSE","Cash.SD","dep.var")
o40.sfl <- mutate(o40.sfl, age.group = "Over 40", tmt = "SFL")

u40.pd <- cbind(u40.pd, u40.sds, pps.out)
colnames(u40.pd) <- c("ATE","RSE","Cash.SD","dep.var")
u40.pd <- mutate(u40.pd, age.group = "Under 40", tmt = "PD")

o40.pd <- cbind(o40.pd, o40.sds, pps.out)
colnames(o40.pd) <- c("ATE","RSE","Cash.SD","dep.var")
o40.pd <- mutate(o40.pd, age.group = "Over 40", tmt = "PD")

pps.het <- rbind(u40.sfl, o40.sfl, u40.pd, o40.pd)
pps.het <- mutate(pps.het, ATE.std = ATE/Cash.SD, RSE.std = RSE/Cash.SD, 
                  CI.high = ATE.std + 1.96*RSE.std, CI.low = ATE.std - 1.96*RSE.std)
pps.het$dep.var[which(pps.het$dep.var=="PPS Index")] <- " PPS Index"

sfl.pps.het.plot <- ggplot(data = pps.het[which(pps.het$tmt=="SFL"),], 
                           aes(y = dep.var, x = ATE.std, group = age.group,
                               fill = age.group, color = age.group, 
                               xmin = CI.low, xmax = CI.high)) + 
  geom_point(position = position_dodge(0.5)) + xlim(-0.5, 0.75) +
  geom_errorbar(width = 0.1, position = position_dodge(0.5)) + 
  geom_vline(xintercept = 0, color = "black") + 
  labs(x = "Treatment Effect (Standard Deviations)", 
       y = "Outcome", title = "Storytelling") + 
  scale_fill_brewer(palette="Paired", direction = -1, aesthetics = c("colour","fill")) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), 
                     axis.text.x = element_text(angle = 45, hjust = 1)) 

pd.pps.het.plot <- ggplot(data = pps.het[which(pps.het$tmt=="PD"),], 
                          aes(y = dep.var, x = ATE.std, group = age.group,
                              fill = age.group, color = age.group, 
                              xmin = CI.low, xmax = CI.high)) + 
  geom_point(position = position_dodge(0.5)) + xlim(-0.5, 0.75) +
  geom_errorbar(width = 0.1, position = position_dodge(0.5)) + 
  geom_vline(xintercept = 0, color = "black") + 
  labs(x = "Treatment Effect (Standard Deviations)", 
       y = "Outcome", title = "Soft Skills") + 
  scale_fill_manual(values=c("#C11B17","#F778A1"), aesthetics = c("colour","fill")) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), 
                     axis.text.x = element_text(angle = 45, hjust = 1))

ggsave(filename="Figures/PPS_Het_SFL.eps", plot=sfl.pps.het.plot, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)

ggsave(filename="Figures/PPS_Het_PD.eps", plot=pd.pps.het.plot, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)

#####Figure OA9 - Heterogeneity by age, speaking up#####
u40.sfl <- as.data.frame(rbind(c(su.het$coefficients[2], su.het$rse[2])))
o40.sfl <- as.data.frame(rbind(c(su.het.2$coefficients[2], su.het.2$rse[2])))

u40.pd <- as.data.frame(rbind(c(su.het$coefficients[3], su.het$rse[3])))
o40.pd <- as.data.frame(rbind(c(su.het.2$coefficients[3], su.het.2$rse[3])))

u40.sds <- as.data.frame(rbind(c(sqrt(mean(cash.u40$fu.su,na.rm=T)*(1-mean(cash.u40$fu.su,na.rm=T))))))
o40.sds <- as.data.frame(rbind(c(sqrt(mean(cash.o40$fu.su,na.rm=T)*(1-mean(cash.o40$fu.su,na.rm=T))))))

su.out <- as.data.frame(c("Speaking Up"))

u40.sfl <- cbind(u40.sfl, u40.sds, su.out)
colnames(u40.sfl) <- c("ATE","RSE","Cash.SD","dep.var")
u40.sfl <- mutate(u40.sfl, age.group = "Under 40", tmt="SFL")

o40.sfl <- cbind(o40.sfl, o40.sds, su.out)
colnames(o40.sfl) <- c("ATE","RSE","Cash.SD","dep.var")
o40.sfl <- mutate(o40.sfl, age.group = "Over 40", tmt = "SFL")

u40.pd <- cbind(u40.pd, u40.sds, su.out)
colnames(u40.pd) <- c("ATE","RSE","Cash.SD","dep.var")
u40.pd <- mutate(u40.pd, age.group = "Under 40", tmt = "PD")

o40.pd <- cbind(o40.pd, o40.sds, su.out)
colnames(o40.pd) <- c("ATE","RSE","Cash.SD","dep.var")
o40.pd <- mutate(o40.pd, age.group = "Over 40", tmt = "PD")

su.het <- rbind(u40.sfl, o40.sfl, u40.pd, o40.pd)
su.het <- mutate(su.het, ATE.std = ATE/Cash.SD, RSE.std = RSE/Cash.SD, 
                 CI.high = ATE.std + 1.96*RSE.std, CI.low = ATE.std - 1.96*RSE.std,
                 grp = paste(age.group, tmt),
                 Treatment = if_else(tmt=="PD","Soft Skills","Storytelling"))

su.het.plot <- ggplot(data = su.het, aes(y = ATE.std, x = Treatment, group = age.group,
                                         fill = age.group, color = age.group, 
                                         ymin = CI.low, ymax = CI.high)) + 
  geom_point(position = position_dodge(0.5)) + 
  geom_errorbar(width = 0.1, position = position_dodge(0.5)) + 
  geom_hline(yintercept = 0, color = "black") + 
  labs(x = "Treatment", 
       y = "Treatment Effect \n (Std deviations)") + 
  scale_fill_brewer(palette="Paired", direction = -1, aesthetics = c("colour","fill")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), 
                     axis.text.x = element_text(angle = 45, hjust = 1)) 

ggsave(filename="Figures/SU_Het.eps", plot=su.het.plot, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)

#####Figure OA10 - Heterogeneity by age, SPS#####
u40.sfl <- as.data.frame(rbind(c(pss.het$coefficients[2], pss.het$rse[2]),
                               c(fss.het$coefficients[2], fss.het$rse[2]),
                               c(ss.het$coefficients[2], ss.het$rse[2]),
                               c(v.het$coefficients[2], v.het$rse[2]),
                               c(r.het$coefficients[2], r.het$rse[2]),
                               c(lid.het$coefficients[2], lid.het$rse[2]),
                               c(con.het$coefficients[2], con.het$rse[2]),
                               c(rmt.het$coefficients[2], rmt.het$rse[2]),
                               c(rm.het$coefficients[2], rm.het$rse[2]),
                               c(goal.het$coefficients[2], goal.het$rse[2]),
                               c(rma.het$coefficients[2], rma.het$rse[2]),
                               c(gc.het$coefficients[2], gc.het$rse[2]),
                               c(sps.het$coefficients[2], sps.het$rse[2])))

o40.sfl <- as.data.frame(rbind(c(pss.het2$coefficients[2], pss.het2$rse[2]),
                               c(fss.het2$coefficients[2], fss.het2$rse[2]),
                               c(ss.het2$coefficients[2], ss.het2$rse[2]),
                               c(v.het2$coefficients[2], v.het2$rse[2]),
                               c(r.het2$coefficients[2], r.het2$rse[2]),
                               c(lid.het2$coefficients[2], lid.het2$rse[2]),
                               c(con.het2$coefficients[2], con.het2$rse[2]),
                               c(rmt.het2$coefficients[2], rmt.het2$rse[2]),
                               c(rm.het2$coefficients[2], rm.het2$rse[2]),
                               c(goal.het2$coefficients[2], goal.het2$rse[2]),
                               c(rma.het2$coefficients[2], rma.het2$rse[2]),
                               c(gc.het2$coefficients[2], gc.het2$rse[2]),
                               c(sps.het2$coefficients[2], sps.het2$rse[2])))

u40.pd <- as.data.frame(rbind(c(pss.het$coefficients[3], pss.het$rse[3]),
                              c(fss.het$coefficients[3], fss.het$rse[3]),
                              c(ss.het$coefficients[3], ss.het$rse[3]),
                              c(v.het$coefficients[3], v.het$rse[3]),
                              c(r.het$coefficients[3], r.het$rse[3]),
                              c(lid.het$coefficients[3], lid.het$rse[3]),
                              c(con.het$coefficients[3], con.het$rse[3]),
                              c(rmt.het$coefficients[3], rmt.het$rse[3]),
                              c(rm.het$coefficients[3], rm.het$rse[3]),
                              c(goal.het$coefficients[3], goal.het$rse[3]),
                              c(rma.het$coefficients[3], rma.het$rse[3]),
                              c(gc.het$coefficients[3], gc.het$rse[3]),
                              c(sps.het$coefficients[3], sps.het$rse[3])))

o40.pd <- as.data.frame(rbind(c(pss.het2$coefficients[3], pss.het2$rse[3]),
                              c(fss.het2$coefficients[3], fss.het2$rse[3]),
                              c(ss.het2$coefficients[3], ss.het2$rse[3]),
                              c(v.het2$coefficients[3], v.het2$rse[3]),
                              c(r.het2$coefficients[3], r.het2$rse[3]),
                              c(lid.het2$coefficients[3], lid.het2$rse[3]),
                              c(con.het2$coefficients[3], con.het2$rse[3]),
                              c(rmt.het2$coefficients[3], rmt.het2$rse[3]),
                              c(rm.het2$coefficients[3], rm.het2$rse[3]),
                              c(goal.het2$coefficients[3], goal.het2$rse[3]),
                              c(rma.het2$coefficients[3], rma.het2$rse[3]),
                              c(gc.het2$coefficients[3], gc.het2$rse[3]),
                              c(sps.het2$coefficients[3], sps.het2$rse[3])))

u40.sds <- as.data.frame(rbind(c(sd(cash.u40$fup.ss,na.rm=T)),
                               c(sd(cash.u40$fuf.ss,na.rm=T)),
                               c(sd(cash.u40$fu.ss,na.rm=T)),
                               c(sd(cash.u40$fu.vew,na.rm=T)),
                               c(sd(cash.u40$fu.rew,na.rm=T)),
                               c(sd(cash.u40$fu.lead,na.rm=T)),
                               c(sd(cash.u40$fu.tg, na.rm=T)),
                               c(sd(cash.u40$fu.rmtime,na.rm=T)),
                               c(sqrt((mean(cash.u40$fu.rm,na.rm=T))*(1-(mean(cash.u40$fu.rm,na.rm=T))))),
                               c(sqrt((mean(cash.u40$fu.goal,na.rm=T))*(1-(mean(cash.u40$fu.goal,na.rm=T))))),
                               c(sqrt((mean(cash.u40$fu.rmasp,na.rm=T))*(1-(mean(cash.u40$fu.rmasp,na.rm=T))))),
                               c(sqrt((mean(cash.u40$com.goal,na.rm=T))*(1-(mean(cash.u40$com.goal,na.rm=T))))),
                               c(sd(cash.u40$fu.sps.ind, na.rm=T))))

o40.sds <- as.data.frame(rbind(c(sd(cash.o40$fup.ss,na.rm=T)),
                               c(sd(cash.o40$fuf.ss,na.rm=T)),
                               c(sd(cash.o40$fu.ss,na.rm=T)),
                               c(sd(cash.o40$fu.vew,na.rm=T)),
                               c(sd(cash.o40$fu.rew,na.rm=T)),
                               c(sd(cash.o40$fu.lead,na.rm=T)),
                               c(sd(cash.o40$fu.tg, na.rm=T)),
                               c(sd(cash.o40$fu.rmtime,na.rm=T)),
                               c(sqrt(mean(cash.o40$fu.rm,na.rm=T)*(1-mean(cash.o40$fu.rm,na.rm=T)))),
                               c(sqrt(mean(cash.o40$fu.goal,na.rm=T)*(1-mean(cash.o40$fu.goal,na.rm=T)))),
                               c(sqrt(mean(cash.o40$fu.rmasp,na.rm=T)*(1-mean(cash.o40$fu.rmasp,na.rm=T)))),
                               c(sqrt(mean(cash.o40$com.goal,na.rm=T)*(1-mean(cash.o40$com.goal,na.rm=T)))),
                               c(sd(cash.o40$fu.sps.ind, na.rm=T))))

sps.out <- as.data.frame(c("SSS - Past","SSS - Future",
                           "SSS - Current","Self-Value Score",
                           "Resilience Score","ID as Leader Score",
                           "Connectedness to Workshop Group","Time to be Role Model","Has Role Model",
                           "Has a Goal", "Aspires to be Role Model",
                           "Achieved Last Goal","SPS Index"))

u40.sfl <- cbind(u40.sfl, u40.sds, sps.out)
colnames(u40.sfl) <- c("ATE","RSE","Cash.SD","dep.var")
u40.sfl <- mutate(u40.sfl, age.group = "Under 40", tmt="SFL")

u40.pd <- cbind(u40.pd, u40.sds, sps.out)
colnames(u40.pd) <- c("ATE","RSE","Cash.SD","dep.var")
u40.pd <- mutate(u40.pd, age.group = "Under 40", tmt="PD")

o40.sfl <- cbind(o40.sfl, o40.sds, sps.out)
colnames(o40.sfl) <- c("ATE","RSE","Cash.SD","dep.var")
o40.sfl <- mutate(o40.sfl, age.group = "Over 40", tmt="SFL")

o40.pd <- cbind(o40.pd, o40.sds, sps.out)
colnames(o40.pd) <- c("ATE","RSE","Cash.SD","dep.var")
o40.pd <- mutate(o40.pd, age.group = "Over 40", tmt="PD")

sps.het <- rbind(u40.sfl, o40.sfl, u40.pd, o40.pd)
sps.het <- mutate(sps.het, ATE.std = ATE/Cash.SD, RSE.std = RSE/Cash.SD, 
                  CI.high = ATE.std + 1.96*RSE.std, CI.low = ATE.std - 1.96*RSE.std)

sps.het$dep.var <- factor(sps.het$dep.var, 
                          levels = c("SPS Index","Achieved Last Goal","Has a Goal",
                                     "Time to be Role Model","Aspires to be Role Model",
                                     "Has Role Model","ID as Leader Score*","Resilience Score", 
                                     "Self-Value Score", "Connectedness to Workshop Group", 
                                     "SSS - Future","SSS - Past","SSS - Current"))

sfl.sps.het.plot <- ggplot(data = sps.het[which(sps.het$tmt=="SFL" & !is.na(sps.het$dep.var)),], 
                           aes(y = dep.var, x = ATE.std, group = age.group,
                               fill = age.group, color = age.group, 
                               xmin = CI.low, xmax = CI.high)) + 
  geom_point(position = position_dodge(0.5))  +
  geom_errorbar(width = 0.1, position = position_dodge(0.5)) + 
  geom_vline(xintercept = 0, color = "black") + 
  labs(x = "Treatment Effect (Standard Deviations)", 
       y = "Outcome", title = "Storytelling") + 
  scale_fill_brewer(palette="Paired", direction = -1, aesthetics = c("colour","fill")) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), 
                     axis.text.x = element_text(angle = 45, hjust = 1)) 

pd.sps.het.plot <- ggplot(data = sps.het[which(sps.het$tmt=="PD" & !is.na(sps.het$dep.var)),], 
                          aes(y = dep.var, x = ATE.std, group = age.group,
                              fill = age.group, color = age.group, 
                              xmin = CI.low, xmax = CI.high)) + 
  geom_point(position = position_dodge(0.5)) +
  geom_errorbar(width = 0.1, position = position_dodge(0.5)) + 
  geom_vline(xintercept = 0, color = "black") + 
  labs(x = "Treatment Effect (Standard Deviations)", 
       y = "Outcome", title = "Soft Skills") + 
  scale_fill_manual(values=c("#C11B17","#F778A1"), aesthetics = c("colour","fill")) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), 
                     axis.text.x = element_text(angle = 45, hjust = 1)) 

ggsave(filename="Figures/SPS_Het_SFL.eps", plot=sfl.sps.het.plot, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)

ggsave(filename="Figures/SPS_Het_PD.eps", plot=pd.sps.het.plot, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)

#####Figure OA11 - Heterogeneity by age, economic outcomes#####
u40.sfl <- as.data.frame(rbind(c(mue.het$coefficients[2], mue.het$rse[2]),
                               c(c.het$coefficients[2], c.het$rse[2]),
                               c(i.het$coefficients[2], i.het$rse[2]),
                               c(ai.het$coefficients[2], ai.het$rse[2]),
                               c(e.het$coefficients[2], e.het$rse[2])))

o40.sfl <- as.data.frame(rbind(c(mue.het2$coefficients[2], mue.het2$rse[2]),
                               c(c.het2$coefficients[2], c.het2$rse[2]),
                               c(i.het2$coefficients[2], i.het2$rse[2]),
                               c(ai.het2$coefficients[2], ai.het2$rse[2]),
                               c(e.het2$coefficients[2], e.het2$rse[2])))

u40.pd <- as.data.frame(rbind(c(mue.het$coefficients[3], mue.het$rse[3]),
                              c(c.het$coefficients[3], c.het$rse[3]),
                              c(i.het$coefficients[3], i.het$rse[3]),
                              c(ai.het$coefficients[3], ai.het$rse[3]),
                              c(e.het$coefficients[3], e.het$rse[3])))

o40.pd <- as.data.frame(rbind(c(mue.het2$coefficients[3], mue.het2$rse[3]),
                              c(c.het2$coefficients[3], c.het2$rse[3]),
                              c(i.het2$coefficients[3], i.het2$rse[3]),
                              c(ai.het2$coefficients[3], ai.het2$rse[3]),
                              c(e.het2$coefficients[3], e.het2$rse[3])))

u40.sds <- as.data.frame(rbind(c(sd(cash.u40$fu.mue,na.rm=T)),
                               c(sd(cash.u40$fu.log.exp.pc,na.rm=T)),
                               c(sd(cash.u40$fu.income,na.rm=T)),
                               c(sqrt(mean(cash.u40$fu.work_lastmonth,na.rm=T)*(1-mean(cash.u40$fu.work_lastmonth,na.rm=T)))),
                               c(sd(cash.u40$fu.e.ind, na.rm=T))))

o40.sds <- as.data.frame(rbind(c(sd(cash.o40$fu.mue,na.rm=T)),
                               c(sd(cash.o40$fu.log.exp.pc,na.rm=T)),
                               c(sd(cash.o40$fu.income,na.rm=T)),
                               c(sqrt(mean(cash.o40$fu.work_lastmonth,na.rm=T)*(1-mean(cash.o40$fu.work_lastmonth,na.rm=T)))),
                               c(sd(cash.o40$fu.e.ind, na.rm=T))))

e.out <- as.data.frame(c("Marginal Utility of Expenditure", "Log(Cons. Exp. per capita)",
                         "Cash Earnings","Earn Any Cash","Econ Index"))

u40.sfl <- cbind(u40.sfl, u40.sds, e.out)
colnames(u40.sfl) <- c("ATE","RSE","Cash.SD","dep.var")
u40.sfl <- mutate(u40.sfl, age.group = "Under 40", tmt="SFL")

u40.pd <- cbind(u40.pd, u40.sds, e.out)
colnames(u40.pd) <- c("ATE","RSE","Cash.SD","dep.var")
u40.pd <- mutate(u40.pd, age.group = "Under 40", tmt="PD")

o40.sfl <- cbind(o40.sfl, o40.sds, e.out)
colnames(o40.sfl) <- c("ATE","RSE","Cash.SD","dep.var")
o40.sfl <- mutate(o40.sfl, age.group = "Over 40", tmt="SFL")

o40.pd <- cbind(o40.pd, o40.sds, e.out)
colnames(o40.pd) <- c("ATE","RSE","Cash.SD","dep.var")
o40.pd <- mutate(o40.pd, age.group = "Over 40", tmt="PD")

e.het <- rbind(u40.sfl, o40.sfl, u40.pd, o40.pd)
e.het <- mutate(e.het, ATE.std = ATE/Cash.SD, RSE.std = RSE/Cash.SD, 
                CI.high = ATE.std + 1.96*RSE.std, CI.low = ATE.std - 1.96*RSE.std)
e.het$dep.var[which(e.het$dep.var=="Econ Index")] <- " Econ Index"

sfl.e.het.plot <- ggplot(data = e.het[which(e.het$tmt=="SFL"),], 
                         aes(y = dep.var, x = ATE.std, group = age.group,
                             fill = age.group, color = age.group, 
                             xmin = CI.low, xmax = CI.high)) + 
  geom_point(position = position_dodge(0.5)) + xlim(-0.65, 0.8) + 
  geom_errorbar(width = 0.1, position = position_dodge(0.5)) + 
  geom_vline(xintercept = 0, color = "black") + 
  labs(x = "Treatment Effect (Standardd Deviations)", 
       y = "Outcome", title = "Storytelling") + 
  scale_fill_brewer(palette="Paired", direction = -1, aesthetics = c("colour","fill")) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), 
                     axis.text.x = element_text(angle = 45, hjust = 1)) 

pd.e.het.plot <- ggplot(data = e.het[which(e.het$tmt=="PD"),], 
                        aes(y = dep.var, x = ATE.std, group = age.group,
                            fill = age.group, color = age.group, 
                            xmin = CI.low, xmax = CI.high)) + 
  geom_point(position = position_dodge(0.5)) + xlim(-0.65, 0.8) + 
  geom_errorbar(width = 0.1, position = position_dodge(0.5)) + 
  geom_vline(xintercept = 0, color = "black") + 
  labs(x = "Treatment Effect (Standard Deviations)", 
       y = "Outcome", title = "Soft Skills") + 
  scale_fill_manual(values=c("#C11B17","#F778A1"), aesthetics = c("colour","fill")) + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5), 
                     axis.text.x = element_text(angle = 45, hjust = 1)) 

ggsave(filename="Figures/Econ_Het_SFL.eps", plot=sfl.e.het.plot, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)

ggsave(filename="Figures/Econ_Het_PD.eps", plot=pd.e.het.plot, device="eps",
       height=2.75, width=6.5, units="in", dpi=400)

#####Table OA5 - Disaggregated Communication Index#####
table <- matrix(NA, nrow=11, ncol=7)

table[1,] <- c("","(1)","(2)","(3)","(4)","(5)","(6)")
table[2,] <- unlist(c("Storytelling", round(creg1$coefficients[2],2),round(creg2$coefficients[2],2),round(creg3$coefficients[2],2),
                      round(creg4$coefficients[2],2),round(creg5$coefficients[2],2),round(creg6$coefficients[2],2)))
table[3,] <- c("",paste0("(",round(creg1$rse[2],2),")"),paste0("(",round(creg2$rse[2],2),")"),paste0("(",round(creg3$rse[2],2),")"),
               paste0("(",round(creg4$rse[2],2),")"),paste0("(",round(creg5$rse[2],2),")"),paste0("(",round(creg6$rse[2],2),")"))
table[4,] <- unlist(c("Soft Skills", round(creg1$coefficients[3],2),round(creg2$coefficients[3],2),round(creg3$coefficients[3],2),
                      round(creg4$coefficients[3],2),round(creg5$coefficients[3],2),round(creg6$coefficients[3],2)))
table[5,] <- c("",paste0("(",round(creg1$rse[3],2),")"),paste0("(",round(creg2$rse[3],2),")"),paste0("(",round(creg3$rse[3],2),")"),
               paste0("(",round(creg4$rse[3],2),")"),paste0("(",round(creg5$rse[3],2),")"),paste0("(",round(creg6$rse[3],2),")"))
table[6,] <- unlist(c("Baseline Level", round(creg1$coefficients[4],2),round(creg2$coefficients[4],2),round(creg3$coefficients[4],2),
                      round(creg4$coefficients[4],2),round(creg5$coefficients[4],2),round(creg6$coefficients[4],2)))
table[7,] <- c("",paste0("(",round(creg1$rse[4],2),")"),paste0("(",round(creg2$rse[4],2),")"),paste0("(",round(creg3$rse[4],2),")"),
               paste0("(",round(creg4$rse[4],2),")"),paste0("(",round(creg5$rse[4],2),")"),paste0("(",round(creg6$rse[4],2),")"))
table[8,] <- unlist(c("N",creg1$N, creg2$N, creg3$N, creg4$N, creg5$N, creg6$N))
table[9,] <- c("Adj. R2",round(as.numeric(summary(creg1)[15]),2), round(as.numeric(summary(creg2)[15]),2),
               round(as.numeric(summary(creg3)[15]),2),round(as.numeric(summary(creg4)[15]),2), round(as.numeric(summary(creg5)[15]),2),
               round(as.numeric(summary(creg6)[15]),2))
table[10,] <- unlist(c("S=SS (p-value)", round(linearHypothesis(creg1, "SFL = PD")[2,4],2), 
                       round(linearHypothesis(creg2, "SFL = PD")[2,4],2),
                       round(linearHypothesis(creg3, "SFL = PD")[2,4],2),round(linearHypothesis(creg4, "SFL = PD")[2,4],2), 
                       round(linearHypothesis(creg5, "SFL = PD")[2,4],2),
                       round(linearHypothesis(creg6, "SFL = PD")[2,4],2)))

myxtable <- xtable(table)

align(myxtable) <- c("l","p{3.55cm}",rep("p{1.6cm}",6))

caption(myxtable) <- "Average Treatment Effects for Individual Outcomes in the Communication Index"

label(myxtable) <- "IndCom"

print(myxtable, floating = TRUE, caption.placement="top",sanitize.text.function = identity,
      include.colnames=F,include.rownames=F,table.placement="H",
      hline.after=NULL,
      add.to.row=list(
        pos = list(0,1,7,nrow(table)),
        command = c(
          paste0("\\toprule"),
          " \\midrule ",
          " \\midrule ",
          "\\bottomrule \\multicolumn{7}{l}{\\multirow{2}{16cm}{\\small Notes: Average treatment effects on individual variables that make up our effective communication index. Columns (1), (2), (4), (5), and (6) show effects on questions that respondents answer on a 5-point frequency scale. Column (1) shows effects on the statement ``I make suggestions to others''. Column (2) shows effects on the statement ``if I have a problem or a new idea that would affect or benefit my community, I raise it to others.'' Column (3) shows effects on the likelihood of speaking up in a group situation in the last two weeks where the respondent did not know everyone present. Column (4) shows effects on the statement ``people understand me when I make suggestions for how to accomplish a task.'' Column (5) shows effects on the statement ``people understand what I am saying when I ask them to do something.'' Column (6) shows effects on the statment ``people understand what I am saying when I ask themto do something.'' We report robust standard errors in parentheses. The bottom row shows the p-values associated with hypothesis tests of equality between the estimated average effect for Storytelling (S) versus Soft Skills (SS).}} \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ "
        )),
      type = "latex",file="IndCom.tex")


